{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.optimize import fsolve\n",
    "from scipy.optimize import root\n",
    "from scipy.optimize import least_squares\n",
    "import math\n",
    "from scipy.misc import derivative\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "rc('text',usetex=True)\n",
    "from pandas import DataFrame\n",
    "import matplotlib.cm as cm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Baseline model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "tolerance = 1e-8 \n",
    "\n",
    "class Baseline:\n",
    "    def __init__(self,\n",
    "                 given_beta = 0.9894,      # discount factor \n",
    "                 given_alpha_b = 0.1,      # acceptance of both monies \n",
    "                 given_alpha_c = 0.65, \n",
    "                 given_c_ub = 0.15,\n",
    "                 given_c_lb = 0.0001,\n",
    "                 given_s = 0.54,           # chance of having alpha_cH when state is H\n",
    "                 given_eta = 0.7,          # prob of getting signals \n",
    "                 given_chi = 0.07,         # prob of getting public signals \n",
    "                 given_B = 0.0098,         # DM utility function parameter \n",
    "                 given_i = 0.01,         # targeted inflation of govt money, (1-beta*(1+i)>0\n",
    "                ):\n",
    "        self.beta = given_beta\n",
    "        self.alpha_b = given_alpha_b \n",
    "        self.alpha_c = given_alpha_c\n",
    "        self.c_ub = given_c_ub\n",
    "        self.c_lb = given_c_lb\n",
    "        self.s = given_s \n",
    "        self.eta = given_eta \n",
    "        self.chi = given_chi \n",
    "        self.B = given_B \n",
    "        self.i = given_i \n",
    "\n",
    "    def equm_allocation(self):\n",
    "        beta = self.beta\n",
    "        alpha_b = self.alpha_b \n",
    "        alpha_c = self.alpha_c\n",
    "        c_lb = self.c_lb\n",
    "        c_ub = self.c_ub\n",
    "        s = self.s\n",
    "        eta = self.eta\n",
    "        chi = self.chi\n",
    "        B = self.B\n",
    "        i = self.i\n",
    "        \n",
    "        zeta = eta - chi\n",
    "        ell_lowerbar = 1e-9\n",
    "                \n",
    "        def utility(q): \n",
    "            return B*math.log(max(q,ell_lowerbar)) - B*math.log(ell_lowerbar) + q\n",
    "        \n",
    "        def S(ell): # q = ell \n",
    "            return B*math.log(max(ell,ell_lowerbar)) - B*math.log(ell_lowerbar)\n",
    "        \n",
    "        def S_prime(ell):\n",
    "            return B/ell\n",
    "\n",
    "        q_star = math.inf \n",
    "        w_star = q_star # when buyer has all the bargaining power\n",
    "\n",
    "        \n",
    "        # closed form solution for pi_d\n",
    "        c_H = chi*(s*c_lb + (1-s)*c_ub)\n",
    "        c_L = chi*(s*c_ub + (1-s)*c_lb)\n",
    "        \n",
    "        if c_L<c_H:\n",
    "            print(\"Error: Wrong c_H and c_L, c_L<c_H\") \n",
    "            \n",
    "        if B*alpha_c - (1+i)*c_H < 0:\n",
    "            print(\"Error: Wrong c_H and c_L ,B*alpha_c - (1+i)*c_H < 0\") \n",
    "\n",
    "        # Do not require phi_hat(0) to be positive or negative    \n",
    "#        if beta*(B*alpha_c - (1+i)*c_L)/(1-beta*(1+i)) > c_L:\n",
    "#            print(\"Error: wrong c_L\",c_L)\n",
    "        \n",
    "        Gamma = ( (1-beta*(1+i)*(1-eta))/(2*s*beta*(1+i)*eta)\n",
    "                - math.sqrt(((1-beta*(1+i)*(1-eta))/(2*s*beta*(1+i)*eta))**2 - (1-s)/s))\n",
    "\n",
    "        k = (math.log(Gamma) / math.log((1-s)/s))\n",
    "        print('k',k)\n",
    "        \n",
    "        if k<= 0: \n",
    "            print('Error: k<=0')\n",
    "            \n",
    "        pi_d_cf = ((c_L*(1-beta*(1+i))/(beta*eta)+(B*alpha_c-(1+i)*c_L)*(Gamma*s + s-1))/\n",
    "            ((B*alpha_c-(1+i)*c_L)*(2*s-1)-(c_H-c_L)*((1+i)*s*(1-Gamma)+(1-beta*(1+i))/(beta*eta))))\n",
    "        \n",
    "        if pi_d_cf >=1:\n",
    "            print(\"Error: pi_d is larger than 1\")\n",
    "            \n",
    "        if pi_d_cf <0:\n",
    "            print(\"Error: pi_d is negative\")\n",
    "        \n",
    "        def update(y): \n",
    "            \" belief updating after receiving a good news, get pi_prime given pi \"\n",
    "            return y*s/(y*s + (1-y)*(1-s))       \n",
    "        \n",
    "        diff = 10 #to calculate the time horizon\n",
    "        pi = pi_d_cf \n",
    "        count = 2\n",
    "        while diff > tolerance:\n",
    "            pi_prime = update(pi)\n",
    "            count += 1 \n",
    "            diff = 1 - pi_prime\n",
    "            pi = pi_prime\n",
    "        n = count #time horizon grid number\n",
    "        \n",
    "        a = np.zeros(n)\n",
    "        pi_array = np.copy(a) \n",
    "        pi_array[1] = pi_d_cf\n",
    "        sbar_array = np.copy(a) \n",
    "        c_array = np.copy(a) \n",
    "            \n",
    "        pi = np.copy(pi_d_cf)  \n",
    "        sbar_array[1] = pi_d_cf*s + (1-pi_d_cf)*(1-s) \n",
    "        sbar_array[0] = 1-s\n",
    "        c_array[0] = c_L\n",
    "        c_array[1] = pi_d_cf*c_H + (1-pi_d_cf)*c_L\n",
    "        \n",
    "        for j in range(1,n-1):\n",
    "            pi_prime = update(pi) \n",
    "            pi_array[j+1] = pi_prime \n",
    "            sbar_array[j+1] = pi_prime*s + (1-pi_prime)*(1-s)  \n",
    "            c_array[j+1] = pi_prime*c_H + (1-pi_prime)*c_L\n",
    "            pi = pi_prime   \n",
    "        \n",
    "        phi_h = (beta*B*alpha_c - beta*(1+i)*c_H)/(1-beta*(1+i))\n",
    "            \n",
    "        phi_l = (beta*B*alpha_c - beta*(1+i)*c_L)/(1-beta*(1+i))\n",
    "                \n",
    "        phi_hat_array = beta*(B*alpha_c - (1+i)*c_array)/(1-beta*(1+i)) # stationary price \n",
    "        \n",
    "        phi_hat_d = beta*(B*alpha_c - (1+i)*c_array[1])/(1-beta*(1+i))\n",
    "        \n",
    "        def eq_ell_b(ell):\n",
    "            return alpha_b*S_prime(ell) - i\n",
    "        \n",
    "        ell_b = fsolve(eq_ell_b, 10)\n",
    "            \n",
    "        # solve for continuous phi \n",
    "        n_c = 3001\n",
    "        pi_c_array = np.zeros(n_c)\n",
    "        pi_c_array[1:] = np.linspace(pi_d_cf, 1, num=n_c-1) \n",
    "        pi_c_array[0] = 0\n",
    "        sbar_c_array = np.zeros(n_c)\n",
    "        c_c_array = np.zeros(n_c)\n",
    "        \n",
    "        for j in range(n_c):\n",
    "            pi_j = pi_c_array[j]\n",
    "            sbar_c_array[j] = pi_j*s + (1-pi_j)*(1-s)  \n",
    "            c_c_array[j] = pi_j*c_H + (1-pi_j)*c_L\n",
    "        \n",
    "        phi_c_hat_array = beta*(B*alpha_c - (1+i)*c_c_array)/(1-beta*(1+i)) # continuous stationary price \n",
    "        \n",
    "        lambda_c_array = np.zeros(n_c)\n",
    "        phi_c_array_cf = np.zeros(n_c)\n",
    "        phi_c_array_cf[:3] = 0\n",
    "        w_c_array_cf = np.zeros(n_c) \n",
    "        w_c_array_cf[:3] = alpha_b*S(ell_b)+alpha_c*S(0)\n",
    "\n",
    "        \n",
    "        for j in range(2,n_c): # j = 1 = d\n",
    "            pi_j = pi_c_array[j]\n",
    "            phi_c_hat_j = phi_c_hat_array[j]\n",
    "            lambda_c_array[j] = (- ((1- beta*(1+i))/(beta*eta)) * phi_hat_d * (pi_j/pi_d_cf)\n",
    "                    *((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k )\n",
    "            #phi_c_array_cf[j] = phi_c_hat_j + (beta*eta/(1-beta*(1+i)))*lambda_c_array[j]\n",
    "            phi_c_array_cf[j] = phi_c_hat_j - phi_hat_d*(pi_j/pi_d_cf)*((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k\n",
    "            w_c_array_cf[j] = alpha_b*S(ell_b)+alpha_c*S((phi_c_array_cf[j]-beta*B*alpha_c)/(beta*(1+i)))-c_c_array[j]\n",
    "                                          \n",
    "        sed_derivative_c = np.zeros(n_c)\n",
    "        ft_derivative_c = np.zeros(n_c)\n",
    "        for j in range(0,n_c):\n",
    "            ft_derivative_c[j] = -(c_H - c_L)/(1-beta*(1+i)) +phi_hat_d*pi_c_array[j]/((beta*(1+i))*pi_d_cf)*((1-pi_c_array[j])*pi_d_cf/(pi_c_array[j]*(1-pi_d_cf)))**k*(k/(1-pi_c_array[j]) -1)\n",
    "            sed_derivative_c[j] = -(phi_hat_d*pi_c_array[j]/(beta*(1+i)*pi_d_cf)\n",
    "                                        *((1-pi_c_array[j])*pi_d_cf/(pi_c_array[j]*(1-pi_d_cf)))**k\n",
    "                                        *((k-(1-pi_c_array[j]))**2-k*pi_c_array[j])/((1-pi_c_array[j])**2*pi_c_array[j])\n",
    "                                   )\n",
    "            \n",
    "        # closed form solution for phi      \n",
    "        \n",
    "        lambda_array = np.zeros(n)\n",
    "        phi_array_cf = np.zeros(n)\n",
    "        \n",
    "        for j in range(2,n): # j = 1 = d\n",
    "            pi_j = pi_array[j]\n",
    "            phi_hat_j = phi_hat_array[j]\n",
    "            lambda_array[j] = (- ((1- beta*(1+i))/(beta*eta)) * phi_hat_d * (pi_j/pi_d_cf)\n",
    "                    *((1-pi_j)*pi_d_cf/(pi_j*(1-pi_d_cf)))**k )\n",
    "            phi_array_cf[j] = phi_hat_j + (beta*eta/(1-beta*(1+i)))*lambda_array[j]\n",
    "        \n",
    "        phi_array_cf[phi_array_cf < 0] = 0 \n",
    "        \n",
    "        # check incentive constraints  \n",
    "        \n",
    "        if not c_array[1]+tolerance >= eta*sbar_array[1]*phi_array_cf[2]:\n",
    "            print('Error: incentive constraint is wrong')\n",
    "        \n",
    "        d = 1 \n",
    "        \n",
    "        # solve for w\n",
    "        \n",
    "        ell_c = np.zeros(n)\n",
    "        w_array = np.zeros(n)\n",
    "        ell_c[-1] = phi_h - c_H\n",
    "        w_array[-1] = alpha_b*S(ell_b) + alpha_c*S(ell_c[-1]) - c_H\n",
    "        sed_derivative = np.zeros(n)\n",
    "\n",
    "        for j in range(0,n-1):\n",
    "            if j <= d: \n",
    "                w_array[j] = alpha_b*S(ell_b) + alpha_c*S(0)  # no need to pay c\n",
    "            else:    \n",
    "                ell_c[j] = (eta*sbar_array[j]*phi_array_cf[j+1] \n",
    "                        +eta*(1-sbar_array[j])*phi_array_cf[j-1]\n",
    "                        +(1-eta)*phi_array_cf[j] - c_array[j])\n",
    "                w_array[j] = alpha_b*S(ell_b) + alpha_c*S(ell_c[j]) -c_array[j]\n",
    "                sed_derivative[j] = (w_array[j] - w_array[j-1])/(pi_array[j] - pi_array[j-1])\n",
    "        \n",
    "        avew = np.zeros(n) \n",
    "        for j in range(n):\n",
    "            avew[j] = pi_array[j]*w_array[-1] + (1-pi_array[j])*w_array[0]\n",
    "            \n",
    "        # solve for closed form welfare \n",
    "        # cannot have closed form solution for continuous welfare\n",
    "        \n",
    "        m1 = (1-beta*(1-eta))/(2*s*beta*eta) - (((1-beta*(1-eta))/(2*s*beta*eta))**2 - (1-s)/s)**0.5\n",
    "        m2 = (1-beta*(1-eta))/(2*s*beta*eta) + (((1-beta*(1-eta))/(2*s*beta*eta))**2 - (1-s)/s)**0.5\n",
    "        \n",
    "        if not m2>1>m1>0:\n",
    "            print('Error with m1 and m2')\n",
    "        \n",
    "        n1 = m1*s/(1-s)\n",
    "        n2 = m2*s/(1-s)\n",
    "        \n",
    "        if not n2>n1:\n",
    "            print('Error with n1 and n2') \n",
    "                    \n",
    "        Omega_cf_array = np.zeros(n) \n",
    "        Omega_d = w_array[0]/(1-beta)\n",
    "        Omega_cf_array[0:d+1] = Omega_d\n",
    "        Omega_cf_array[-1] = w_array[-1]/(1-beta)\n",
    "\n",
    "#         for j in range(d+1, n):\n",
    "#             sum_bw = 0\n",
    "#             pi_j = pi_array[j]\n",
    "#             for i in range(d, 1000):\n",
    "#                 if i<n:\n",
    "#                     pi_i = pi_array[i]\n",
    "#                     w_i = w_array[i]\n",
    "#                 else: \n",
    "#                     pi_i = pi_array[-1]\n",
    "#                     w_i = w_array[-1]\n",
    "                \n",
    "#                 if i <= j:\n",
    "#                     b_ji = ((pi_j*(m1**(j-d))/(beta*eta*(1-s)*(m2-m1))) *(m2/(m1**(i-d-1))-m1/(m2**(i-d-1)))/pi_i)\n",
    "#                     sum_bw = sum_bw + b_ji*w_i\n",
    "#                 else: \n",
    "#                     b_ji = pi_j*m1*((m2**(j-d))-(m1**(j-d)))/(beta*eta*(1-s)*(m2-m1)*(m2**(i-d-1))*pi_i)\n",
    "#                     sum_bw = sum_bw + b_ji*w_i\n",
    "#             Omega_cf_array[j] = Omega_d + sum_bw\n",
    "        \n",
    "        w_c_array = np.zeros(n_c)\n",
    "        ell_c_c = np.zeros(n_c)\n",
    "        ell_c_c[-1] = phi_h - c_H\n",
    "        w_c_array[-1] = alpha_b*S(ell_b) + alpha_c*S(ell_c_c[-1]) - c_H\n",
    "\n",
    "        for j in range(0,n_c-1):\n",
    "            if j <= d: \n",
    "                w_array[j] = alpha_b*S(ell_b) + alpha_c*S(0)  # no need to pay c\n",
    "            else:    \n",
    "                ell_c_c[j] = (eta*sbar_c_array[j]*phi_c_array_cf[j+1] \n",
    "                        +eta*(1-sbar_c_array[j])*phi_c_array_cf[j-1]\n",
    "                        +(1-eta)*phi_c_array_cf[j] - c_c_array[j])\n",
    "                w_c_array[j] = alpha_b*S(ell_b) + alpha_c*S(ell_c_c[j]) -c_c_array[j]\n",
    "\n",
    "        avew_c = np.zeros(n) \n",
    "        for j in range(n):\n",
    "            avew_c[j] = pi_c_array[j]*w_c_array[-1] + (1-pi_c_array[j])*w_c_array[0]\n",
    "            \n",
    "        # iteration to solve for welfare    \n",
    "        Omega = np.zeros(n) \n",
    "        Omega[0:d+1] = Omega_d \n",
    "        Omega[-1] = w_array[-1]/(1-beta)\n",
    "        \n",
    "        Omega_test = np.copy(Omega)\n",
    "        Omega_array = np.copy(Omega)\n",
    "        diff = 10 \n",
    "        count = 1\n",
    "        while diff > tolerance and count <= 3000: \n",
    "            \n",
    "            for j in range(d+1, n-1): \n",
    "                sbar_j = sbar_array[j]\n",
    "                w_j = w_array[j]\n",
    "                #Omega_j = min(max((w_j + eta*sbar_j*beta*Omega[j+1] \n",
    "                #               + eta*(1-sbar_j)*beta*Omega[j-1] \n",
    "                #            + (1-eta)*beta*Omega[j]),Omega[0]),Omega[-1])\n",
    "                Omega_test[j] = max((w_j + eta*sbar_j*beta*Omega[j+1] \n",
    "                               + eta*(1-sbar_j)*beta*Omega[j-1] \n",
    "                            + (1-eta)*beta*Omega[j]),Omega[0])\n",
    "                 \n",
    "            count += 1 \n",
    "            if count ==3000: \n",
    "                print(\"Count=3000 number of iterations reached at welfare \")\n",
    "            diff = np.max(np.abs(Omega - Omega_test))\n",
    "            Omega = np.copy(Omega_test)\n",
    "            \n",
    "        Omega_array = np.copy(Omega) # update welfare array \n",
    "        \n",
    "        aveOmega = np.zeros(n)\n",
    "        \n",
    "        # full disclosure, pi W_H + (1-pi) W_L\n",
    "        for j in range(n): \n",
    "            pi_j = pi_array[j]\n",
    "            aveOmega[j] = pi_j*Omega_array[-1] + (1-pi_j)*Omega_d\n",
    "                          \n",
    "        return phi_c_hat_array,ell_c_c,pi_c_array, phi_c_array_cf,w_c_array_cf, pi_array, pi_d_cf, phi_array_cf, w_array, Omega_cf_array,Omega_array,aveOmega, avew, lambda_c_array, phi_hat_array\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k 1.2247234465321786\n",
      "k"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/liangfan/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/minpack.py:175: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
      "  improvement from the last ten iterations.\n",
      "  warnings.warn(msg, RuntimeWarning)\n",
      "/var/folders/68/f_1qhmls3_b2w5ws_5ln_kdh0000gn/T/ipykernel_48523/3293611824.py:167: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  ft_derivative_c[j] = -(c_H - c_L)/(1-beta*(1+i)) +phi_hat_d*pi_c_array[j]/((beta*(1+i))*pi_d_cf)*((1-pi_c_array[j])*pi_d_cf/(pi_c_array[j]*(1-pi_d_cf)))**k*(k/(1-pi_c_array[j]) -1)\n",
      "/var/folders/68/f_1qhmls3_b2w5ws_5ln_kdh0000gn/T/ipykernel_48523/3293611824.py:167: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ft_derivative_c[j] = -(c_H - c_L)/(1-beta*(1+i)) +phi_hat_d*pi_c_array[j]/((beta*(1+i))*pi_d_cf)*((1-pi_c_array[j])*pi_d_cf/(pi_c_array[j]*(1-pi_d_cf)))**k*(k/(1-pi_c_array[j]) -1)\n",
      "/var/folders/68/f_1qhmls3_b2w5ws_5ln_kdh0000gn/T/ipykernel_48523/3293611824.py:169: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  *((1-pi_c_array[j])*pi_d_cf/(pi_c_array[j]*(1-pi_d_cf)))**k\n",
      "/var/folders/68/f_1qhmls3_b2w5ws_5ln_kdh0000gn/T/ipykernel_48523/3293611824.py:168: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  sed_derivative_c[j] = -(phi_hat_d*pi_c_array[j]/(beta*(1+i)*pi_d_cf)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 1.0646560321906267\n",
      "k 1.0646560321906267\n",
      "k 1.0528353962349624\n"
     ]
    }
   ],
   "source": [
    "eta_grid = [0.2,0.7,0.8,0.99]\n",
    "\n",
    "for j in range(len(eta_grid)): \n",
    "    eta = eta_grid[j]\n",
    "    result = Baseline(given_eta = eta).equm_allocation() \n",
    "\n",
    "    globals()[f\"phi_c_hat_array_{j}\"],globals()[f\"ell_c_c_{j}\"],globals()[f\"pi_c_array_{j}\"],globals()[f\"phi_c_array_cf_{j}\"],globals()[f\"w_c_array_cf_{j}\"],globals()[f\"pi_array_{j}\"], globals()[f\"pi_d_cf_{j}\"], globals()[f\"phi_array_cf_{j}\"],globals()[f\"w_array_{j}\"],globals()[f\"Omega_cf_array_{j}\"],globals()[f\"Omega_array_{j}\"], globals()[f\"aveOmega_{j}\"], globals()[f\"avew_{j}\"],globals()[f\"lambda_c_array_{j}\"],globals()[f\"phi_hat_array_{j}\"]= result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABILElEQVR4nO3deViU1dvA8e9hdXcEFRNFHdxwSUW0tFJTsMwlMzRbzDbxZ/tiku2LZlhWlllivdmiaS5lLqlgpqa5AOaWO66oKcIICsh23j9mQERAUZgZZu7PdXk5zzZzzwPMPec857mP0lojhBBC2BsXWwcghBBCFEUSlBBCCLskCUoIIYRdkgQlhBDCLkmCEkIIYZfcbB1AWaldu7Zu3LixrcMQQghRSrGxsYla6zqF1ztMgmrcuDExMTG2DkMIIUQpKaUOF7VeuviEEELYJUlQQggh7JIkKCGEEHZJEpQQQgi7JAlKCCFEqZlSUpj/4ZNsXjSt3F7DYUbxCSGEKF8XLlzg21k/k3P2CP2Tv6PJyeOkqP7Qf2S5vJ5TJaizZ8+SmJhIZmamrUMR18jDw4PatWtTs2ZNW4cihNO4cOEC8UkXmDFrFh+OHsWEXp6492qJ9xNzCOzUt9xe12kSVEZGBv/99x8NGjSgcuXKKKVsHZIoJa016enpHDt2DE9PTypVqmTrkIRwaFprunS7nRTtyZgeNRjvGk3fsDq0HvI63j2eBNfyTSFOk6BOnz5NnTp1qFKliq1DEddIKUWVKlWoXbs2p0+fpmHDhrYOSQiHs2HDBqKjo+k//Gm+WLmbKm5ZDK32Dw+7KbI6PEq3kNehipdVYnGaBJWRkUG9evVsHYYoA9WrV+fMmTO2DkMIh3Hs2DHq16+Pi4sLP/6ylOlTPmXjGfiw5s+07H6YbL9bcek7EU+f1laNy2kSVHZ2Nm5uTvN2HZqbmxvZ2dm2DkMIh/DHH38QHBzMlB/m83dGfXZkNuLH529jsPuH5FavD3d+h1uru8EGl0Wcapi5XHdyDPJzFOLaZWRk8MILLzB79mwAXOs1p3W/x4lYl0zbY7PZ4P0moZ4b4dYXcHkmBloPtElyAiu1oJRSBsBo+ddJax1exD6hgAkwaq0ji1snhBCidEwmE/v27aNTp054enry559/koYHS843Ye2+RHp1asfEyl/ifW4v+N0Od30ItZvZOmyrtaCGAEFa63kASqmwghstiQitdbRlObiodVaKVQghHMrw4cMZNGgQOTk5/HsihQ7PTGW5ZzdOJBxhRZOf+CbnNbxdzsPg72DYL3aRnMBKCUprHVmgBWQE4gvt0qnAunggsJh1l1BKhSmlYpRSMadPny77wEWxwsPD8ff3x9/fn8jIkhu38fHxhISEUKtWravaXwhxfZYtW0bbtm0xmUwAvPnmm3z2zY8889M/9P3sL7YcSeb7NluIcn+R5v8tg1tfgKc327Q7ryhWHTWglDICSXmtogIMhZa9i1l3CUvSiwQICgrSZROluJLw8HCio6OJjY0FoGPHjnh5eREaGlrk/iEhIUybNo3g4GCio6MZPHgwAGFhYUXuL4QonczMTH777Tc6duxIkyZNqF27NnXq1OHUqVOczXFn5n5Xfv3nHJXd0/kg6DyDT03Gdf8OMNpPd15RrD1IIlRrXVRNDBNQeGB9UeuEHYiMjGTu3LkYDAYMBgPh4eFMmDChyH3j4uIwGo0EB5t7aIODgxk7diwRERHWDFkIh5SRkQFAUlISQ4cO5fvvvwcgKCiIHxYs4dsdGfSctJqlO07w3M01iWu7gKE7RuCakWx33XlFsVoLSikVqrWeaHkcqLWOK7B5MxdbTEYgyrJceF2ZemfRTv49nlLWT1sqrerX4K3+1r234HrExcVhMpkwGo3564xGI3FxcUXuHxgYyLRplxaTNBqNxMcX7uUVQpTGXXfdRfXq1ZkzZw716tVj48aNtG/fnlMpGUz98wCzNh4BYPhNvrxgWEv19RMhK83cndftZfCoauN3cGXWGsUXDEQopcZaVoVb1kdprUO01vOUUmMs+xkKDIy4bJ0wi4yMJDY2loiICKKjzacmKiqKkSNHEhh42eW6MpOUlITBYLhkXV6yMplMl20ruD3P5s2byzVGIRzRpk2bWLJkCe+88w5g7o0oWO7LP6AtE1fs5bv1h8jK0QwJasCLLZOps/pJ2GL/3XlFsUqCsiQX/yLWhxR4PNHyMLqkdWWpIrVcCoqOjmbIkCHExsYyePBgoqIuNi4nTJjA3LlzizwuPDw8/6JpcQwGA2PHji0y0QBFHu/lZe6JLSp5FXV8XhehEKJkCQkJ+Pj44Obmxrp165g8eTJPPvkkPj4+vPjiiwCkZ+bw7fqDfPXnAVIvZDOwvS8vdDHgFxsBP8+CGr7m7jwb3Wx7PaS0QgVkNBoxGAzExMQwffr0/PXx8fH5yaIoZXHdx2AwFJvkSnrtPL169WL69On516SEEEVbt24d3bp1Y+HChfTr14+wsDBGjhyZX080KyeXuTHHmLxyL/+lXKBXy7q83NuflkfnwszxFa47ryiSoCqgvC6z+Pj4S7rK8rr4ylNRSSjvetKVWk8hISGMHDmy2NF+QjizzMxMXnnlFdq1a8fw4cPp3Lkzb7/9NjfeeCMAVauak4zWmqXbTzJpxR7iE8/TsVEtpjwQSCe3ePjtbji5rUJ25xVFElQFFR0dTVBQ0GXrSuo6K4suvsDAQAwGA3FxcfnJMTo6+ootopCQEAYPHixDy4Uo4OzZs/z777906dIFDw8P1q9fT+XKlQFwd3fnjTfeuGT/dfsTiVi2m23HztLcpxrTHw4iuLEHauU7EDsDqteDwTOg1cAK151XFElQFVRUVBQhIfmX8PITRkmtmLIa2h0WFkZ4eDhRUVHEx8czYcKES7oaTSbTJa27wYMHS3ISoghhYWGsWrWKhIQE3N3dWbduHa6urpftt/3YWSYu383afYnUr1mJD0NvZFAHX1y3z4Ypb0B6Etw8CnqMhUo1bPBOyodTFYt1JPHx8Zd0lc2ZM8dq13UiIiIIDAzE39+fkJAQIiIiLoklMjIy/2bcuLg45s2bx8iRI1FKXfJPCGezcuVKWrVqRV7lm1dffZUlS5bkz7RQODkdTDzPU7Pi6D/lL3YknOX1vgH8MboHg/3O4fp9P/h1FHg1gbDVcOcEh0pOAEprxyjAEBQUpGNiYordvmvXLgICAqwYkXV17NiRiIgIpxl84Og/T+EYMjMzWbRoEW3btqV58+bs3LmTp59+mqlTp5b4+3sqJYPJK/cxZ/NR3F1deOK2JozoZqSGSyasngh/TwGPahDyDnR4GFwqdltDKRWrtQ4qvF66+BxEXFyc0yQnIexdRkYGlSpVIjU1lQceeIAXX3yRCRMm0Lp1a1atWlXscakZWUxbHc83fx0kKyeX+zv78UyvptStXgl2L4Hfw+HsUWj/kDk5Va1txXdlfZKgHEBcXJyMjBPCTgwcOBCtNQsXLsTb25sNGzbkj8QrTlZOLj9tOsLk6H2cOZ9J/3b1eSmkOY1rV4XkwzArHPb+DnVbwaPLoFEXK70b25IE5QACAwPlxlchbCQ2Npb58+czfvx4lFKEhISQm5ubv71Dhw7FHqu1ZvnO/5i4bDfxiee5qYkX3/YN4MYGBsjOhLWTYPWHoFwg5D3zQAhXdyu8K/sgCUoIIUopISGB2rVr4+npyebNm/niiy8YNWoUDRs25Kmnnrqq59hyJJn3l+5i86Fk/OtU5euHg+gVUNc8gOjgWljyEiTugZb9oE8E1GxQzu/K/lTsK2tCCGFlMTEx+Pn5sWjRIsA8GeCJEydo2LDhVR1/5EwaT82K456p6zmYmMb4e9qw/PluBLfyQZ0/DQvC4Lt+kJ0BD/wMQ2c6ZXICaUEJIUSJcnJyGDNmDC1atCAsLIwOHTrw7rvv0qlTJ4D8G2uvJPl8Jp//sZ8fNhzCzcWFZ3s1I6ybkWqebpCbA5u/hpXvQmYa3DYabnsJPKqU51uze5KghBCikLNnz7J9+3ZuvfVWXF1diY2NxcUylNvV1ZXXXnvtqp8rIyuH79YfYsqq/Zy/kM2QoIa8ENIcnxqWSuTHt8DiF+F4HDTpBn0/rvAlisqKJCghhMA8YCHvBvJnn32WhQsXcvLkSSpVqsTKlSuLrPBQktxczW9bj/Ph8j0kmNK5vUUdXukTQIt61c07ZJyFP8aZW05VasOgr6FtqEOUKCorcg1KCOH0Vq9eTatWrUhISABgzJgxrFixAk9PT+DyCg9Xsv5AInd/sY7n5/yDoYo7M5+4iW8f7WxOTlrDtrnweZA5OXV6Ap7eDDcOluRUiLSghBBOJysri0WLFtGiRQtat25NgwYN8PHxISkpCV9fX1q3vra54vafOseEpbtYufsUvobKfHJfO+5u54uLiyXxJO4zj847uBrqd4AHfzb/L4okLShxTcLDw/H398ff35/IyMjr3r+0zyfEtUhPT8//f9iwYXzzzTcA+Pv78+eff9K2bdtret7k85m8/dtO7vx0DZsOJhF+Z0tWvtSdezo0MCenrHRzd96XXeH4P9B3EjyxUpLTFUgLSpRaeHg40dHRxMbGAuY6gF5eXsVWswgPDyc+Pp4DBw5gMplo0qQJRqMxvzTTlbYLURbuu+8+kpOTWbFiBTVq1GD9+vXX3FLKk5WTyw9/H2byyn2kZmRxf2c/XghpTu1qnhd32hcNS14E02G4cSj0fg+q1b3Od+MktNYO8a9jx466JP/++2+J28XVMxgM+sCBA/nL06ZN04GBgUXum5ycrIFL9h8zZowODg6+qu3FkZ+nuJK4uDg9evRonZubq7XWOjIyUn/66af5y9cjNzdXR+08qW//cJVuFL5YP/T1Br37RMqlO6Wc0Prn4Vq/VUPrzzpqHb/6ul/XUQExuojPdat18SmlQpVSUcVsC1RKHVBKxVr+RVjWJyulopRSY6wVpyhZXFwcJpMpf1ZfMM/wGxcXV+T+eRXmC+7fqVOn/PVX2i5EaRw/fjy/G2/btm1MmzYtf8bnESNG8Nxzz133VC+7TqTw0DcbeeL7GFDwf48E8f1jnS+OzsvNgU3TYUon2L0Ubn8NRq0zDyEXpWK1Lj6t9TylVHHzkXtprf3BnKwAk2X9YK11dLkF9fsrcHJ7uT39VanXFvp8UOrDIiMjiY2NJSIiguho8ynKm/K94DTwZS0pKemySRHzkovJZLritO95+5c0s++VtgtRlO3bt9O+fXtmzJjBsGHDGDp0KKGhoflTpV+v06kX+DhqD3M2H6VGZXfe7t+KB29uhLtrge/5J7bB4uchIRaMPcz3NHn7l8nrOyO7uAZVKAkZtdbzLI8NSimj1jq+qOOUUmFAGICfn185R2k/oqOjGTJkCLGxsQwePJioqIsN0wkTJhRbOLYspnwv6ngvLy+g6OSVNy19weSV1zoymUxX3H41CU84p9zcXMLDw/Hz8+OZZ56hTZs2jB8/nltvvRUAT0/P/GHi1yMjK4f/W3eQqasOkJGVwyNdm/Bsr6YYqnhc3OnCOfhzAmz4Eqp4yT1NZcQuElQepVSY1rrgEC4vIEkpNU1rfVnry7JvJJgnLCz1C15Dy8UeGI1GDAYDMTExl0y1Hh8fn58silIWU74bDIZik1xRr20wGAgODs5PpPHx8flx5CWfK20XIk9KSgpbtmyhe/fuuLi4sGPHDjIzMwFQSvHKK6+U2WtprVm6/SQTft/FseR0ggPq8updARjrVLt0x91LYOkYSDkGHR+F4Legcq0yi8OZ2dsw85CCC1rrSK21CTAppWTCI4u8LrX4+PhLuvOioqIICQkp7rAyUVQSyuvjLy6h5LXo/P39iYiIIDw8/JK4r7RdODddYNbvMWPG0LdvX86dOwfAkiVLmDx5cpm/5tajJoZM+5unZsVRzdONmU/cxNfDO12anM4eg58egNkPmKdaf2wF9P9UklMZspsWlFLKUGg5DPPIjqKvvju56Ojo/O6xgutKmheqLLr4AgMDMRgMxMXF5SeR6OjoEoeEGwyGS7ohBw8ezMiRI696u3Be69ev57HHHmPp0qUYjUZefPFFHn300fzrSi5lPNX5ibPpfLhsDwu2JFC7mgcTBrVlSFBDXF0KdNXlZMPGr2DV+6BzIfgd6PKUU83TZC1WS1BKqWAgSCkVmneNSSkVpbXO+8rvBSQVOORnwJjXcipwXUpweWspL2GU1C1WFl18AGFhYYSHh+d3yU2YMOGSrkaTyXRJ6y4uLi6/W3LixInEx8cTFhZ2SewlbRfOIysri8WLF9OkSRPat2+Pn58f9erV4+zZswA0b968XF43IyuH6WvimfrnAXK0ZlQPf57s4U/1SoWSzrFYWPyceXBVszvgrg+hVqNyiUlYsYtPax2tta5VMNEUSE5oreMLXmfSWpu01nFa63la63BrxVlRxMfHX3Jj7Jw5c6x2Y2tERASBgYH4+/sTEhJCRETEJbFERkYyePDg/OXo6GiaNGlCrVq1OHDgQP4Nvle7XTi+tLQ0wJygHnnkkfxqIg0aNODPP/8scVba66G1ZtmOEwR/vJpJUXvp0aIOK1/sTvidLS9NThlnzSWKvu4F5xNhyPfwwBxJTuVMFezfrciCgoJ0SffO7Nq1i4CAACtGZF0dO3YkIiLCaaovOPrP05kMGzaMQ4cOsXbtWgB27NhBy5YtcXMr3w6ePSdTeWfRTtYfOEMLn+q81b8VXZvWvnQnrWHnAlg2Fs6fhs5h5vuaKtUo19icjVIqVmsdVHi93VyDEtcnLi7OaZKTqNi2bt3Kt99+y6RJk3B1dSU4OJgzZ87kT3fRpk2bcn19U1omn0Tt5ceNR6jm6cY7A1rz4E1+uLkW6lBKOghLR8P+aLihHdw/G3xl8I41SYJyAHFxccXWwRPCHpw4cYIaNWpQtWpV9uzZw//93/8xcuRIAgICGD58uFViyMnVzNp0hI9X7OFsehYP3OTHSyEtqFXV49IdszNh/Wew5kNwcYM7P4BOI8BVPi6tTc64AwgMDCxx9J4QtrR3715atWrFl19+yYgRIxg4cCB9+/YtswoPV2ND/Bne/m0nu0+mcrPRi7f6tybghiK66Q7/ba4EcXo3BPSHOyOgpq/V4hSXkgQlhChTWmvGjh2Lt7c3L7/8Ms2aNeODDz6gV69eAHh4eODh4XGFZykbx5LTmLB0N0u2n8DXUJmpDwbSp029y+vxpSVB1Juw5Qeo2RDunwMt7rRKjKJ4kqCEENctNTWVzZs307NnT5RS7NmzhxtuuAEwV3gYPXq0VeNJz8zhq9UH+Gr1AZSCF4KbM7K7kUruhWbG1Rq2zoYVr0G6Cbo+Az3Ggof1WneieJKghBDXJG9QA8Abb7zBV199xYkTJ6hVqxbz588v85torzamJdtPMGHpbhJM6fS78QbG3hWAr6Hy5Tsn7jPP03RwDTToBP0+MRdvFnbD3kodCSEqgE2bNtGyZUt27doFwDPPPMOqVavybxS3RXL693gKQyM38PSsLdSo7M7ssJuZ8kDg5ckpKwNWTbDMbrvVXHH8sRWSnOyQtKCEEFeUnZ3N4sWL8fX1pVOnTjRu3Jj69evn18Tz9/fH398200qcTc/i4xV7+GHDYWpWdmfcwDbc39nv0vJEeeJXw+IXIOkAtAmFO96H6j7WD1pcFUlQQohinT9/nqpVq5Kbm8uIESMYMGAAnTp1om7duqxatcqmseXmaubFHiNi2W6S0zJ56OZGvBTSgppViqiJdz4Rlr8G22ZDrSbw0AJo2sv6QYtSkQQlhCjS448/zrZt29i8eTMeHh6sXr263Grhldb2Y2d587cdbDliomOjWnx/d2da1695+Y5awz8zYcXr5jmbbhsN3UaDexHXpITdkQQlhADMJYamTZvGxx9/jLu7OyEhIQQEBJCTk4OrqyutWrWydYiY0jL5cPkeZm06gndVDz4a3I5BHXxxKao7L3G/+Z6mQ2uh4c3mqTDqSnmsikQSlIOYOHEigYGBUu5IlMrJkyepUqUKNWrU4ODBg8yYMYPHH3+c9u3bM3ToUFuHly83VzMn5igTl+0mJSObR7o25vng5tSsXER3XnYmrJtsrgThVsk8Oi/wEbDBwA1xfeQn5gAiIyMJDg4mPj6e6OhoW4cjKojDhw/ToEEDZsyYAUCfPn04efIk7du3t2lchf1z1MQ9U9cxdsF2mtWtzuJnbuWt/q2LTk6H/4Zpt8GqcdDyLnh6EwQ9JsmpgpIWlAMIDg7GaDQSGBhIXJzM7yiK9/rrr1OpUiVef/11GjVqxKRJk7jrrrsAcHNzK/cK4qWRdD6Tict2MyfmKLWrefLpfe25u339y6tAgPkm2+i3IHaGuRLEAz9D8zusHbIoY/bz2yiuWd4U8IBMlS4ukZqayt9//03v3r0B8zxiVapUyd/+3HPP2Sq0YuUVdf1o+R7OXcjm8Vua8Fxws8snDwTLdBi/wLJXzNNhdHnaXAnCs9rl+4oKRxKUA5k3b55UNReXVHgYP348kyZNIiEhgbp16zJz5syiWyB2IvZwMm/9toMdCSl0MXrzzt2tae5TveidTUfMkwjuW2GeDuOBOVC/fCY2FLYhHbNCOJAtW7YQEBDAli1bABg1ahRr1qyhTp06AHabnBLPXeDluVu598v1nE69wOf3d2DWiJuKTk452bB+CnxxExxaZ77Z9ok/JDk5IGlBOQiTyZRfZkY4j5ycHBYvXkydOnXo2rUrjRs3xtfXl4yMDAAaNWpEo0b2Oy15Xnfeh8t2k5aZw8juRp7t2YyqnsV8NB3fAouegxNbofmdcNeHYPCzbtDCaqzWglJKhSqlokrYnqyUilJKjSl0TLBSKqy84urRo0f+KKasrCx69OjBjz/+CEBaWho9evRgzpw5AJw9e5YePXqwYMECABITE+nRoweLFi0CzEN2e/TowbJlywA4evQoPXr0yB9ZFx8fT48ePVi9ejUAe/bsKbP3ER0dLUPMnUheiSGtNU8++SRffPEFALVq1WLlypV06dLFluFdlR0JZxk0dR1v/LqD1vVrsuz52xjbJ6Do5HThnHna9ek9IfUkDP7OPMOtJCeHZrUWlNZ6nlJqZAm7DNZa54+RVkqFWo6LVkqFKaWCC24XZpGRkQDExsYC5iQ4ZsyYkg4RFdyoUaNYu3Yt27dvx83NjT/++OOSgTL2LiUji0nLzbXzvKpeYXQewJ5l5qnXzx41Dxnv9RZUNlg1ZmEjWmur/QOiStgWChgLLEcAgZbHwcCYIo4JA2KAGD8/P12Sf//9t8TtFU1ycrIODQ3VBw4c0FprPXfuXK211lFRUXrMmDG2DM0qHO3nWZKdO3fqUaNG6fT0dK211vPnz9cTJ07UmZmZNo6sdHJzc/WvW47poHFRuvEri/Ubv27XprQS3kPKCa3nDNP6rRpaT+ms9eG/rRessCogRheRF+xpkIQXkKSUmmZZNhTa7l34AK11pNY6SGsdlHcR2FmMGDGCkSNHYjQaMZlM+d+gjUYj8+bNs3F04nqdPHmS5ORkABISEvjhhx/Yvn07AIMGDeLll1/G3b2IYdd26sDpczz0zUaem/0P9WpUYuFTt/Du3W2Kvtk2Nxc2fwNTOplbTz1fh5Frwe9m6wcubMpuBklorSMBlFImS/eeCXPSEoWYTCaio6OZO3cuYL7+lDe8PD4+vkJ194jLnTx5koYNGzJu3DjCw8Pp1asXJ06coFq1indvT0ZWDl+s2s+01fF4urvw3t2teeCmRkVPhQFwapd5EMTRjdCkG/T7FLxtM42HsD27SFCWQRAxWuuCZRA2c7EVZQSKHWDhbEpKQtOmTWPkyJIu9Ql79M4775Cdnc17771HvXr1mDx5MiEhIYB58r+KmJxW7TnFWwt3ciQpjYHt6/Nq3wDqVq9U9M5Z6bDmI3MNPc/qMPAraDcU7HRYvLAOqyUopVQwEKSUCtVaz7Osi9JahwA/A8YCAyPyto+xHGfQMkAiX3HVIiIjIzEajXKzbgVw7tw51q5dS58+fQBzXbycnJz87U8++aStQrtuJ86m8+6if/l9x0mMdaoy64mb6Nq0dvEHxK82Vx1Piod290Pv8VD1sh594YSU+fpUxRcUFKRjYmKK3b5r1y4CAhyn1H58fDzz5s3DaDQSHx+PwWDAaDQ6zVDzivjzzPtbU0rx9ttv8+6773LkyBEaNGhwSfWHiiorJ5cZ6w7xSfRecnI1z/ZqxhO3NcHTzbXoA86fgRWvwdafwMtorjpu7GHVmIV9UErFaq2DCq+3iy4+UXpGo5ExY8Ywb948GVZeAezYsYPBgwfzzTff0LVrV0aMGEHv3r3x9fUF7LfCw9WKOZTE67/uYPfJVHq2rMs7A1rT0KtK0TtrDVtnw/JX4UIK3PYSdHtZJhEUl5EEJUQ5yM3NZfHixdSsWZPu3bvTuHFjGjRokN+N5+vrm5+cKrKzaVl8sGwXP206yg01K/HVQx25o7VP8Qn3zAFY/AIcXA0NOkP/yeBj+4kQhX2SBCVEGUpNTaV6dXP9uOeff57AwEC6d+9OtWrViIpynHE+Wmt+23qc9xb/S3JaFiNua8Lzwc2LL1GUkwXrP4PVE8HVA/p+DB0flXmaRIkkQVVwMiDCfjz//PMsXryYvXv34uLiwvLly2ncuLGtwypzR5PSeO3XHazZe5obG9Tku8c607p+zeIPOBYLvz0Dp3ZCwADoMxFq3GC9gEWFJQlKiGu0e/duPv30UyZNmkTVqlUJDg6mfv36ZGVl4enpSbNmzWwdYpnKysnlm78O8mn0XlyV4q3+rXi4S+Pi72m6cA7+GAcbv4LqN8DQWdCyr3WDFhWaUyUoRxgpJS6OhrOF//77Dzc3N7y9vUlMTOTHH3/koYce4tZbb6Vfv37069fPZrGVpy1Hkhm7YDu7T6bSu5UP79zdmhtqljCoYV+U+VrT2aPQ6Qlz/bxKNawXsHAITpOg3N3dSU9Pv2Q2UVExpaen26TMz5kzZ/Dz82Ps2LG8/fbb3HLLLZw8ebJC3kR7tVIzsvjQUtjVp7p5EMSdbeoVf8C507B8LGyfC7VbwGPLpUSRuGZOk6Dq1q1LQkICvr6+VK5cWVpSFZDWmvT0dBISEvDx8bHKa44fPx6TycSHH36It7c3n3/+OT169ADMQ8MdOTkt23GSt37bwanUCwzv0piXejcvetp1sAwd/8kydPycedr1W18AN0/rBi0citMkqBo1zN0Lx48fJysry8bRiGvl7u6Oj49P/s+zrJ07d44///wzv6vu+PHjJCUl5XcPh4WV29RkduO4KZ23fttJ1L//EXBDDaYNC6J9Q0PxByQdNFeCiP8TGt4E/T+Dui2tFK1wZE6ToMCcpMrrg01UXAUrPEyZMoWxY8dy4MABjEYjU6ZMcZrWdk6u5rv1h5i0Yg85WjO2T0seu7UJ7q7FDAXPyYYNU2HV++DiBn0nQcfHZOi4KDNOlaCEKGzPnj3cc889fP755/Tq1YtHH32U2267jSZNmgAVv8LD1dqRcJZXf9nOtmNn6d68DuMGtim+EgTA8X9g0bPmqddb9DVPvV6z4t94LOyLJCjhVHJzc1myZAmVK1cmODiYRo0a4efnh4vlW7+Pj4/Vrm/Zg/TMHD6N3svXfx2kVhUPPr+/A/1uvKH4xJyZBn++D39Phaq1Ycj35nubnCSRC+tymmKxwrmlpKRQo0YNtNa0bt0af39/Fi1aZOuwbGr9gUTGLtjO4TNpDO3UkLF9AqhZpYTRkQdWma81JR+CwOEQ8g5UrmWtcIUDk2KxwmmFh4cza9YsDh06hKurK4sWLcLPz8/WYdnM2fQsJizdxezNR2nsXYVZI26iq38J02GkJcHy12DrLPBuCo8sgca3Wi9g4bQkQQmHs2/fPj766CMiIiIwGAwEBwfj5eVFZmYmlStXxt/feWdoXbbjJG8u3MGZ85mM7G7kheDmVHIvZjoMrWH7PFj2CmSY4LbRlqrjxUw6KEQZkwQlHMKpU6cA8/1uKSkpzJw5k/vuu4+ePXsSEhKSPzutszqVmsFbC3fy+46TBNxQg2+Gd6JtgxLq55mOwOIXYX8U+HaEAb+BT2vrBSwEkqCEA0hJSaFx48Y8++yzfPDBBwQGBjp8hYerpbVmbswxxi35l4zsXF6+owVh3YzFDx3PzYFNkbDyPfPynRHQeQS4FNPKEqIcSYISFVJERAQJCQl89tln1KhRgylTptC1a1fA8Ss8XK0jZ9IY+8s21u0/Q+fGXky4ty3+dUo4Lyd3mIeOJ8RCs97m+5oMznutTtieJChRIZw7d47o6GgGDhwImLv0Tpw4kV/h4bHHHrNtgHYkJ1fz7bqDfLRiD24uLowb2IYHOvvhUlzV8awMWDMR1k2GSga49xtoc68MHRc2Z7UEpZQKBUZqrS+7GKCUMgBGy79OWutwy/pkIAaI0lpPtFaswj4UrPDw9ddf88ILL7Bz505atWrFRx995DQ30ZbGrhMpvDJ/G1uPnaVXy7qMu6dNyVXHD66FRc9B0gFo/yD0HgdVvKwXsBAlsFpNEq31vBI2DwGC8vZRSuUVPBustQ6R5OR84uPjadOmDUuWLAFg2LBhrFmzhoCAAMB5KjxcrQvZOUxasYf+n//FseR0Pru/A18PDyo+OaUnmycR/K4f6BwY9isMnCrJSdgVu+ji01pHFlg0AnlzYxuUUkatdbwNwhJWpLVmyZIluLq60qdPHxo2bEiTJk2oVMk8pNnb25vbbrvNxlHapy1Hknl53jb2nzrHoA6+vNGvFbWqehS9s9bw70JY+jKknYFbnoPur4CHTEMj7I9dJKg8SikjkKS1jras8gKSlFLTtNYji9g/DAgDnPrGy4rMZDJhMBhQSvHmm2/i7e1Nnz59cHd3Z/HixbYOz65lZOXwSdRepq+Nx6dGJb59tBO3t6hb/AEpx2HJaNizBG5oBw/NM/8vhJ2yqwQFhBZMRHktK6WUSSkVWrib0LI9Esyljqwaqbhub7zxBpGRkRw9ehQPDw8WLFiAr68UHL0asYeTeXneVuJPn+f+zn68elfL4udqys2FuO8g6k3IyYKQ9+DmJ8HV3v78hbiU3fyGWhLQRMvjQCAIiNFax9k2MlFW4uPjiYiIYNy4cdSpU4fg4GCqVq1KVlYWHh4eNG7c2NYh2r30TPO1pm/WHaR+zcr88HhnbmtWp/gDzhwwD4I4tBYa3wYDPgMvo/UCFuI6WG2QhFIqGAiyjObLWxdVYFuEUipWKRWLuWvvZ8u2ULjiIAthp/KGgwOkpaXx448/EhsbC0D37t155ZVXqFq1qi1DrDA2H0rirs/W8vVfB3mgsx/LX+hWfHLKyTYPG/+yK5zYZp5EcPgiSU6iQpFq5qLcpKWl4ePjw2OPPcbkyZMBOH/+vCSkUkrLzObD5XuYsf4QvobKTLz3Rro2LaG468ntsPBpOPEPtOwHd30ENW6wWrxClFaZVTNXSjUGgrXWX1uW2wPxWuuU6w1SVHyTJk1i7969TJs2jSpVqjB16lSCgi7+3klyKp2N8WcYM38bh8+k8XCXRoTf2ZKqnsX82WZlwJoPYd2n5mkwBn8Hre6WG25FhVWqBKWUegL4H1AT+Nqy2h8YC9xXtqGJiuD8+fMsW7aMQYMGoZQiOTmZU6dOkZubi4uLC8OGDbN1iBXS+QvZTFy2m+/+PoyfVxV+GnEzXfy9iz/gyAbzfU2Je6HdA3DHeLmnSVR4periU0rtwzx4IVZr3bTA+jNa6xL+esqfdPFZj9YarTUuLi5MmzaN//3vf8TGxhIYGJhfekhcu/UHEgmfv41jyek80rUxL9/RgioexXyXvJAKK9+FTdOhZkPo/wk0DbZuwEJcp+K6+Eo7SMJLa322qOe/trBERXP06FHatm3LggULABg6dCirV6+mQ4cOgFR4uB7nL2Tzxq87eGD6RlyVYk5YF97q37r45LQvCr642ZycbhoJT/4tyUk4lNJeg1qplBoE5De7lFJzgOjiDxEV3dKlS8nOzmbAgAHUr1+fpk2b5lcLr1mzJt26dbNxhBXfpoNJjJ67laPJaTx+axNG925BZY9iprg4fwaWj4Vtc6B2C3h8BTTsbN2AhbCC0nbx1QRWAoGYyxEFAfFAL1sPkpAuvrKVnJxMrVq1ALjllltwdXVlzZo1No7K8WRkme9r+vqvgzSsVYVJQ9rRqXEx1460hp0LYOkY8wy3t74I3UaDm6dVYxairJXJKD5L916QUqoX5pp5E7XWK8soRmEnxo0bx8cff0xCQgKVK1fmp59+ol69erYOy+FsO2bixZ+3sv/UOR662Y+xfQKKH6F3NgGWvAR7f4f6gTBgIdRrY92AhbCy0o7iawxgSUorLet6Yh5mfqisgxPWcejQISZMmMCbb76Jr68vvXr1wtXVlezsbEDqHJa1rJxcPv9jP1+s2k+dap58/1hnujUv5obb3FyImwFRb5nLFPUeDzePkhluhVMo7TWoaUAEcKjAulqWdTLMvAI5ffo0Fy5coEGDBmRnZzNz5kz69u2Lr68vXbp0oUuXLrYO0SHtOZnKiz//w87jKQzq4MtbA1pTs3IxNfTOHIDfnoXDf0GTbtB/slSCEE6ltAkqSGv9R8EVWuv5SqnI4g4Q9iczM5MWLVoQGhpKZGQkTZs25dSpU1SpIlMulJecXM30tfF8vGIv1Su58dVDHbmzTTHdpjnZ8PcU+HMCuHrCgM+hwzC54VY4ndImKKWUqq61Ti28vqwCEuVj8uTJbNmyhRkzZuDh4cHUqVNp27Zt/nZJTuXnUOJ5Xpq7ldjDydzR2ofx97SldrViBjac2Aa/PQ0ntkqZIuH0Spug5mLuznsyb4VS6ksshV2F/UhLS2PJkiWEhoailCI1NZWkpCSys7Nxc3Nj6NChtg7R4eXmamZuPMz7S3fj5qr45L52DGzvW/S9YlkZsDrCXOC1ireUKRKC0o/iG2mpOH4G8/Byo+X/XuURnCidghUe5s6dyyOPPMK6devo2rUrr732mtxEa0XHTemMmbeNv/Yn0q15HSLubVv89OuH15vLFJ3ZD+0fhN7jpEyREFxDsVitdUfL9BhNMI/ek2HmduDEiRP07t2bMWPGMGzYMEJDQ2ncuHH+YAdJTtaz8J8EXv91Bzm5mvH3tOGBzn5Fn/8LqRD9Nmz+Ggx+8NACaCrf9YTIc00TFhaYkl3Y0LJlyzh//jz33nsv9erVIyAgAC8v8zfvqlWr0r17dxtH6FzOpmXxxsId/Lb1OB0b1eLjIe1o5F1M9fZ90eaJBFMSzLPb3v4aeFazbsBC2LkSK0kopSYA0/LucbJUMy9S3vQbtuIslSSSkpLyk1BwcDApKSls2rTJxlGJ9QcSeennrZxOvcBzvZoxqoc/bq5FlLpMT4blr8E/M6F2c7j7CylTJJzetVaSGIy5pNEhy/L/itlPc3H6DVFOJk6cyHvvvcfx48epXr06M2bMoG7durYOy6ldyM5h0oq9TF8bTxPvqswf1ZV2DQ1F77x7CSx+Ac4nwm0vQbcx4F7JqvEKUZGUmKAKTqlh0dPWNfecyZEjR3j//fcZM2YMRqORXr16kZ2dTV6rt0GDBjaO0LntOZnKc7O3sPtkKg/e5MdrfQOKrjx+PhF+HwM75oNPG3jgZ6jf3urxClHRlPYa1EGlVOMi7oMSZSQxMZG0tDT8/MwX1mfOnElwcDBGo5GOHTvSsWNHW4fo9HJzNd+uP0TEst3UqOTG/z0SRM+WPpfvmF/c9WXISDFfZ7rleXDzsHrMQlREpU1Q84CJwKjSvpBSKhQYqbUOKWG7CTBqrSOLW+fIsrOzadWqFXfccQc//PADDRs25NSpU1SuXMzwZGF1J89mMHruVv7an0hwQF0+uPfGom+6TT1pLu66e7G5uOvdX4BPK+sHLEQFVtoEtQKYrpQyYr42ZcrbcKVBElrreUqpkUVtsyQitNbRSqkwyzB2Q+F1jjh6cOrUqaxZs4bZs2fj5ubG1KlTadmyZf52SU72Y8m2E7z6y3Yys3OZMKgtQzs1vHz4uNawdTYsewWy0iHkXbj5KXC9pgGzQji10v7VjMV8Y643ULAUwfUOkugEzLE8jsc835R3EesqfIJKS0tj0aJFhIaG4urqSlpaGqmpqVy4cAFPT09CQ0NtHaIoJDUji7d+28mCuATaNTTw6X3taVK7iOHjpqOw+HnYHw0Nb4a7p0DtZlaPVwhHUdpKEpcNAywjhkLL3sWsu4RSKgwIA/ueEkJrTW5uLq6urixevJihQ4eycuVKevbsyUsvvcTo0aNtHaIoRtyRZJ6bvYXjpgye69WMp3s2xb3w8PG8KTFWvAk6B/pMhE4jwKWIYeZCiKt2TX9BSqmeSqknlFK3l1EcJqBwbZei1l1Cax2ptQ7SWgfVqVPMfDo2lpiYSLt27fjmm28AuPvuu1m1ahU9evQApMKDvcrJ1Xyxaj+Dv/obreHnkTfzQkjzy5NT0kH4foB5+LhvB3jyb7hppCQnIcpAaScs7IC5m60Wllp8Sqn9QO/rnLBwMxdbTHnXtwxFrKsQVqxYwZkzZ7j//vvx9vambdu2+PiYR3l5enrmJydhn06ezeD5OVvYEJ9E/3b1GX9PG2pUKjRnU24ObIqEle+CcjXP1RQ4XIq7ClGGSqwkcdnO5mQUp7UeUmDdXKC61vrOKxwbjLka+git9TzLuqi8UX1KqTFAHBCotZ5Y3Lri2LqSxJkzZ/D2NvdC9u/fnyNHjvDPP/9IC6mCWbHzJGPmbyMzO5d3BrQmtGODy3+Gp/eap8Q4uhGa9YZ+n0BNuSdNiGtVXCWJ0iaoJK21V6F1BiBJa23TPg1bJqjJkyczduxYEhISqFWrFsePH6d27dp4eMj9LhVFRlYO45fs4ocNh2njW4PPhnbAWKdQbbycbPj7c1g1AdwrQ58IuPE+aTUJcZ2KS1ClTSrRSqnqhdZpHGB0XWkkJCQwatQodu3aBUDPnj0vmc6ifv36kpwqkD0nU7l7yjp+2HCYEbc1Yf6orpcnp/92wte9zNXHm4XAU5ug3VBJTkKUo9IOM08C/lBKzSmwLgRIVkrlD0XTWn9UFsHZkzNnzpCSkkKTJk1wc3Nj5syZ3HbbbQQEBNC2bdtLZqcVFYPWmh83HGbckl1Ur+TOd491pnvzQoNtsjPhr49hzUdQqSYMngGtBkpiEsIKSpug8ppghadj9Qb8LY814FAJKjc3l3bt2tGlSxfmzp2Lj48P//33n9xEW4Eln89kzPxtRP37H92b1+Gjwe2oU71QRYjjW2Dh0/DfDmg7BO78AKpedreDEKKc2Mt9UHYnMjKS33//nV9++QUXFxemTJmCv79//nZJThXXhvgzPD/7H86cv8DrfQN47JYmuLgUaBEVnH69Wl24fza06GO7gIVwUlJ/xSI9PZ2FCxdy77334u7uTlZWFhcuXOD8+fNUrVqVgQMH2jpEcZ1ycjVTV+3nk+i9NPKuyi/Db6GNb81LdzoWA78+CYl7oMMw8/TrlQ02iVcIZ1eqUXz27HpH8S1atIgBAwawZMkS7rrrrjKMTNiD06kXeGHOP/y1P5G729dn/D1tqeZZ4PtZVgb8+T6s/xyq14cBk6FpsO0CFsKJXOuEhU6jT58+rFq1im7dutk6FFHG1u9P5Lk5/5CSnsUHg9pyX+EirwVbTYEPm1tNlWoW/4RCCKuQBGXh5uYmFR4cTE6u5vM/9jF55T6Mtavyw+OdaVmvxsUdCreaHloATXvZLmAhxCUkQQmHdCo1g+dn/8P6A2cY1MGX9wa2oWrBLj1pNQlh9yRBCYfz175Enp+zhXMXspkYeiODC5YrklaTEBWGJCjhMHJyNZOj9/L5qv3416nGrBE309ynQOETaTUJUaFIghIO4VRKBs/8tIWNB5MI7diAd+9uTRUPy6+3tJqEqJAkQYkKb0P8GZ6etYXzF7L5aHA7QjsWqCx+LAZ+HQWJe83TYfR+T1pNQlQQkqBEhaW1JnJNPBOX76GRVxVmPnETLepZuvSk1SREhScJSlRIKRlZvDx3K8t3/sddbesRce+NVM+bVPCyVtM4qFSj5CcUQtgdSVCiwtl1IoVRP8ZyNDmd1/sG8PitTcyj9KTVJIRDkQQlKpQFccd49Zft1Kjkzk8jbqZzE8v8mdJqEsLhSIISFcKF7BzeXfQvMzce4WajF5/d34G61StJq0kIByYJSti9Y8lpPDUzjq3HzvK/7v6M7t0cN1cXaTUJ4eCslqCUUqGACTBqrSMLbQsE5lq2A0RrrcOVUslADBCltZ5orViF/Vi77zTP/LSFnBzNtGEduaN1Pci+ANETzPM1SatJCIdllQRlSU5oraOVUmFKqWCtdXSBXby01v6WfQO5mKgGF9pPOIm8IeQRy3bT3Kc6Xz7UkSa1q8KJrfDL/+DUv+b5mu54X1pNQjgoFyu9Ticg3vI4HggsuLFQEjJqrfP2NSiljFaIT9iRtMxsnp39DxN+302ftjew4MmuNKnlAX9GwPSekJYED8yFu6dIchLCgVmri89QaNm7qJ2UUmGFuv+8gCSl1DSt9cii9gfCAPz8/MooVGFLR5PSGPF9DHv+SyX8zpb8r7sRdXq3udV04h9oOxj6TIQqXrYOVQhRzqzVgjJhTjZXElJwQWsdqbU2Aaa8bsIitgdprYPq1KlTJoEK2/lrXyL9p/zFcVM6Mx7tzKhujVHrP4Np3eDsURjyPdz7tSQnIZyEtVpQm7nYijICUYV3UEoZCi2HATFa67jyDk7Yltaar9ceZMLvu2hWtzrThnWksToJ3z4IRzdCy37Q71OoJl9ChHAmVmlBaa3nAUalVDBgyLvmpJQqmKi8gKQCyz9b9gkt8BzCwaRn5vDc7H8Yv3QXd7apx4JRN9M4fhZ8dSuc3g33RMJ9P0pyEsIJKa21rWMoE0FBQTomJsbWYYhSOJqUxsgfYtl1MoXRvVvwZHt31G9Pw8E10DQYBnwONerbOkwhRDlTSsVqrYMKr5cbdYVNbDqYxP9+jCUrJ5f/Gx7E7Wkr4MuxgIb+k8033ubNgiuEcEqSoITV/RxzlNd+2U7DWlX4NrQBjdY9DftWQOPb4O4voFYjW4cohLADkqCE1eTkaj74fRfT1x7kVn9vIjscosrs4ZCdAXdGQOcwcLHWwFIhhL2TBCWsIjUji+dm/8Mfu08xKqgGL+dMxWXxQmjQCQZ+CbWb2TpEIYSdkQQlyt2RM2k88f1mDpw+z4wup+ix9znIOAu93oKuz4Kr/BoKIS4nnwyiXG2MP8P/foylSu55/m65kLpbfoF6N8LDC8Gnta3DE0LYMUlQotzM2XyE13/dQb/qB/jQ/UvcDp6EbmOg28vg5mHr8IQQdk4SlChzubmaiOW7mbF6N5NrL6LPuQWoKkZ4fAU0uOxWByGEKJIkKFGmMrJyeOnnrRza8TdrDNPxOXcQgh6H3u+BR1VbhyeEqEAkQYkyk3Q+k5HfbaRTwg98Vmk+Lm614cH50CzY1qEJISogSVCiTBxMPM9r3yzklbRP6Oi+FwIGQr9PpPK4EOKaSYIS1y3m4BmWfh/B1/o7PDzdod9087xNUqpICHEdJEGJ67Ji0zbcFj/Hmy5xpDe4FbfB06BmA1uHJYRwAJKgxDXRWrN8/td02v4O1VwukNZzHFVufUpKFQkhyowkKFFqueln2Tr9f9yZtJQjlZpRbfh3VKkvN90KIcqWJChRKpkH1pD60xPcmHWKv+o/QtfHJuLi7mnrsIQQDkgSlLg62RfIjHoPt41TSMmty9pO/8fA/oNsHZUQwoFJghJXdnoPWT8/hsfpHfyU04vK/SYw8KYWto5KCOHg5Iq2KJ7WsGk6uV914/zpI4zKGUO9B7+S5CSEsAqrtaCUUqGACTBqrSOL2J4MxABRWuuJV3OMKEfnTsPCp2Dfcv5WHXhTjeLDEb0J9Ktl68iEEE7CKi0oS6JBax1tWS6q9s1grXVIoeR0pWNEedgXBV92IffAKt7Xj/Ky++tMG3WXJCchhFVZq4uvExBveRwPBBaxj0EpZSzlMaIsZaXD0jEwM5Rzbl4MyBzPqpr3MP+pW2hat7qtoxNCOBlrJShDoWXvIvbxApKUUtOu9hilVJhSKkYpFXP69OnrDtKpndwOkbfDpmkcbDqcmxJfBZ8A5ozswg01K9s6OiGEE7JWgjJhTkDF0lpHaq1NgKnAtaerOSZIax1Up06dMgrVyeTmwt9fwPSekJ7E2psi6bXzDgIa1GXWiJvxqioTCwohbMNaCWozF1tERiCq4EZLS6hwF16Jx4gykHICfhwEy1+FpiH83GkOw1ZXo6t/bb5/vDM1KrnbOkIhhBOzSoLSWs8DjJaBDoYCAx/yks7PluW8gRHzijtGlJFdi+HLrnB0I/SfzFc3vMuY348THODD18ODqOIht8gJIWxLaa1tHUOZCAoK0jExMbYOw/5lpsHysRA7A25ojx40nY+3aD7/Yz/929Xn4yHtcHeV2+OEENajlIrVWgcVXi9fk53Jf//CvEfh9B645Xn07a8yMfogX/55gPuCGvL+oLa4usgcTkII+yAJyhloDTH/Z77WVKkmDPsFbexBxLI9fLX6AA/e5Md7d7fBRZKTEMKOSIJydOnJ8NszsGsRNA2GgV+hq9aW5CSEsHuSoBzZkQ0w/wlIPQG9x8HNT6GV4oNlu5m2Op6Hbvbj3QGSnIQQ9kkSlCPKzYG1H8OfE8DQEB5fAb4d0VpLchJCVBiSoBxNynFYEAaH1kLbwdD3Y6hUQ5KTEKLCkQTlSPYsg19HQXYG3D0V2j8AypyEPonex7TV8Tx4kyQnIUTFIAnKEWRnQvRbsGEq1GsLod9C7Wb5myPXHOCzlfsY3LGBDIgQQlQYkqAquuTD5nubEmKh80gIeRfcK+Vv/nHDYd5fupu+N97AB/feKMlJCFFhSIKqyHYvhV//BxoY8gO0GnDJ5gVxx3hj4Q56tqzLJ0Pay024QogKRRJURZSTBSvfhfWfwQ3tYPAM8DJessuyHScYPXcrNzfxZuqDgXi4SfkiIUTFIgmqojmbYO7SO7oROj0Bvcdf0qUHsGbvaZ75aQvtGhr4engQldxdbRSsEEJcO0lQFcm+aFgwAnIy4d5voG3oZbtsO2bifz/G4l+nGjMe6UxVT/kRCyEqJvn0qghyss033a79COq2hiHfXTJKL8+hxPM8+u1malXx4PvHOlOzisznJISouCRB2bvU/2D+4+YbbwMfhj4Twf3yKdhPp17g4f/bRK7WfP94Z+rWqFTEkwkhRMUhCcqeHdkIPz8MF1Jg4FfQ/v4idzt3IZtHZ2zidOoFZo24Cf861awcqBBClD1JUPZIa9j8NSx7BWo2hGELwKd1kbtmZucy6sdYdp1IZfrDHengV8vKwQohRPmQBGVvstJh8Quw9SdodgcMioTKhiJ31VrzyoJtrN2XyMR7b6RnSx/rxiqEEOVIEpQ9ST4Ecx6Ckzugx6vQ7WVwKf7+pal/HmBBXALPBzdjSKeG1otTCCGswGoJSikVCpgAo9Y6stA2A2C0/OuktQ63rE8GYoAorfVEa8VqE/ujYd7jgIYH5kDzO0rcfen2E3y4fA8D2tXnuV6Xj+gTQoiKzirlBSzJCa11tGU5uNAuQ4AgrfU8y/Ywy/rBWusQh05OWsOaj+DHUKjhCyNWXTE5bT1q4sWf/yHQz8DE0BtRSkoYCSEcj7Xq33QC4i2P44HAghu11pEFWlXGAvsalFKX1vBxJJlpMO8x+OM9aDMInogCb/8SDzluSueJ72OoXc2TyIelSoQQwnFZK0EZCi17F7WTJRkl5bW0AC8gSSk1rZj9w5RSMUqpmNOnT5dZsFaRchy+7QM7f4Hgt82VITyqlnhIWmY2j38XQ3pmDt8M70Ttap7WiVUIIWzAWgnKhDnZXEmo1npk3oKlZWUCTHndhAVZtgdprYPq1KlTZsGWu2OxEHk7nNkP9/8Et76QP7FgcbTWhM/fzu6TKXz+QAda1KtupWCFEMI2rDVIYjMXW1FGIKrwDkqp0LxrTUqpQCAIiNFax1kpRuvYNhcWPgXVfWBYFPi0uqrDvvnrIIu2HuflO1pwe4u65RykEELYnlVaUJbBD0bL4AhDgcESUZb/g4EIpVSsUioWc2vrZ8u20ALPUXHl5kL0O7DgCWjQCUb8edXJaf2BRCb8vps7W9fjyR4lX6MSQghHobTWto6hTAQFBemYmBhbh1G0rAzzxII7f4HA4XDXR+DmcVWHHjel0+/zv/Cq6sGvT91CNalOLoRwMEqpWK11UOH18mlX3tKSYPYDcORv6D0Oujx9xetNeTKycvjfj7FkZucybVhHSU5CCKcin3jlKfmQ+f4m02EI/dY8lLwUxi35l23HzjJtWEcpACuEcDqSoMpLQhzMGmKenv3hhdCoa6kO/337CX7ccISwbkbuaF2vnIIUQgj7JQmqPOxfaa6pV7U2PLIU6jQv1eFHk9IYM38b7RoaGN27RTkFKYQQ9k0SVFn7d6G5pl6dlvDQfPNw8lLIzM7l6Z+2ADDl/g54uFnrVjUhhLAv8ulXlrbMhLmPgG8gPLK41MkJ4KMVe9h61ETEvTfS0KtK2ccohBAVhCSosrLhK1j4JDTpDsN+KXYOp5Ks3XeayDXxPHSzH3e1vaHsYxRCiApEElRZWPcZLAuHlv3MU2VcoaZeUc6mZfHy3G3416nK632v7gZeIYRwZJKgrteGLyHqDWh9Dwz+DtyurYDr24t2cvrcBT65r71UKBdCCCRBXZ9N02HZKxAwAAZNB9drG3OydPsJftmSwDM9m3JjA0PZxiiEEBWUJKhrFfcDLB0NLe4yT5Xh6n5NT3MqNYPXftlOW9+aPHV70zIOUgghKi5JUNdi73JY9Bz494LBM666rl5hWmte+2UH5zNz+OS+dri7yo9DCCHyyCdiaSXEmoeS12sLQ76/5mtOAMt2nCTq3/94KaQ5TevK/E5CCFGQJKjSSD4MM4dA1Trw4FzwvPb6eGfTs3jrt520uqEGj9/apAyDFEIIxyCVJK5WZhrMeRBys+ChBVDt+iYNnLhsN4nnLvD18CDcpGtPCCEuIwnqamgNi1+AkzvggZ+h9vUNZth8KImZG4/w+K1NZNSeEEIUQ766X43YGbBtNvQYC817X9dTZeXk8uqC7fgaKvNiSOmKyAohhDORFtSVJMXD8tfA2AO6vXzdT/f934fZd+oc0x8OoqpMQCiEEMWSFlRJcnPgl1Hg4gZ3TwWX6ztdiecu8Gn0Xro1r0NwwPVdwxJCCEdnta/wSqlQwAQYtdaRV7P9SseUu7jv4egGGPgV1PS97qf7aPke0jNzeLNfK9RVTvsuhBDOyiotKEuiQWsdbVkOvtL2Kx1T7tKT4Y/3oNEt0G7odT/d9mNnmRNzlOFdG9O0rkzfLoQQV2KtFlQnYI7lcTwQCERfYbv3FY4pU+feqocHmfnLLmhc0Aw8249dr/9+3c+fk6vxquLBs72aXfdzCSGEM7BWgjIUWva+iu1XOgalVBgQZlk8p5Tac23h5asNJF666qnrfMqLDgKGN8vs6cpbEefCKcl5uEjOxUVyLi4qi3PRqKiV1kpQJsCrlNuvdAyW61Jldm1KKRWjtQ4qq+eryORcmMl5uEjOxUVyLi4qz3NhrQS1mYstIiMQdRXbDVc4RgghhAOzyiAJrfU8wGgZ6GAoMPAhqrjtxR0jhBDCOVhtmLnWeqLlYXSBdSFX2H7ZunJm/aHs9kvOhZmch4vkXFwk5+KicjsXSmtdXs8thBBCXDOpJCGEEMIuOWUxuGupauGoSnqvSikD5gEqRqCT1jrc6gFa0dX+3JVSEc5+LpRSgZh/L/KuITss+by4yPJeRxa8PFPEdhNldC6crgV1LVUtrB2jtVzFex0CBOV9AFnuO3NIV/tzt6w3WjE0q7vKczHW8nvhpZRy2PNxFZ8XwUC8ZXu8JXE7rJK+jJTHZ6fTJSjMVSviLY/zKlSUZrsjKfG9aq0jC3wLMhbY1xFd8edu+SB25HOQp8RzYfmislkpZbT8jjjyObnS70UMMDevRam1jrNmcHamzD87nTFBGQotX01VC0dlKLRc5Hu1fDAnOfhQf0Oh5aLOhdHBP4zzGAotFz4X/pZ1SUqpaZauYEdlKLR8ybnQWpuAacBcoKN1QrJbhkLL1/3Z6YwJykTpq1o4KhNX915DtdYjyzkWWzNRwrlQSgU7eIIuyMSVfy8OWD6cY7lYbswRmbjC7wUQrbX2B0x53VxOykQZf3Y6Y4K6lqoWjuqK71UpFZp3P5qD969f6VwkFaiyb3Tyc7G5wGMD5g8mR3WlcxFYoFtvAs7z5bYoZf7Z6XQJ6lqqWtgu2vJ1pXNhWR+hlIpVSsXiwH98V/F7EWdZ58XlXRkO5Sr/Rgx5F8EdeeTalc4FEKmUCrNsH+LI5wLyPxOCCrYUy/OzU27UFUIIYZecrgUlhBCiYpAEJYQQwi5JghJCCGGXJEEJIYSwS5KghBBC2CVJUEIIIeySJCghhBB2SRKUEEIIuyQJSgghhF2SBCWEEMIuSYISwg4ppUKVUgeUUsmF/h1QSkXYOj4hrMEpp3wXwp7lzVCrtfZXSoVprSPzqqc7+YR4wslIC0oI+5NUYGptg+X/YJxjNl8h8kmCEsLOWCYCzG9JWfjnrRfCWUiCEsJ+hQN5c+oYS9pRCEckCUoI+xVc8JpT3gSBQjgLSVBC2CHLjKXzCqyKw8Fn8hWiMJlRVwghhF2SFpQQQgi7JAlKCCGEXZIEJYQQwi5JghJCCGGXJEEJIYSwS5KghBBC2CVJUEIIIeySJCghhBB26f8BQaKRyAceahYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To make figures of asset prices\n",
    "\n",
    "plt.plot(pi_c_array_0,phi_c_array_cf_0,label=r'$\\eta$ = ' + str(eta_grid[0]))\n",
    "#plt.plot(pi_c_array_1,phi_c_array_cf_1,label=r'$\\eta$ = ' + str(eta_grid[1]))\n",
    "#plt.plot(pi_c_array_2,phi_c_array_cf_2,label=r'$\\eta$ = ' + str(eta_grid[2]))\n",
    "plt.plot(pi_c_array_3,phi_c_array_cf_3,label=r'$\\eta$ = ' + str(eta_grid[3]))\n",
    "\n",
    "plt.plot(pi_array_0,phi_hat_array_0,label=r'$\\hat \\phi$',linestyle=':',color='k')\n",
    "\n",
    "#plt.axvline(x=pi_d_cf_0)\n",
    "#plt.axvline(x=pi_d_cf_1)\n",
    "plt.xlabel(r'$\\pi$', fontsize=16)\n",
    "plt.ylabel('price', fontsize=16)\n",
    "plt.ylim([0,phi_c_array_cf_0[-1]*1.03])\n",
    "\n",
    "plt.legend(fontsize=16)\n",
    "plt.tight_layout()\n",
    "\n",
    "# DataOut = np.column_stack((pi_c_array_0,phi_c_array_cf_0))\n",
    "# np.savetxt('phi_loweta_lowcH.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_c_array_3,phi_c_array_cf_3))\n",
    "# np.savetxt('phi_higheta_lowcH.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_c_array_0,phi_c_hat_array_0))\n",
    "# np.savetxt('phi_hat_lowcH.dat', DataOut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3C0lEQVR4nO3dd1yc153v8c+ZGXobEKigPsi2ZEmWDdiWbbkK3dhJ7CQ2oM3azm6SNWyy2U3bFUnWN5vs3qyC0rNOEUpyncQpikh1bpwY3EtsC7AsW5KbRr0iDUMTnXP/ODMIEJ2ZeR5mfu/Xi9dUnvkxQvPllOccpbVGCCGEsBuH1QUIIYQQI5GAEkIIYUsSUEIIIWxJAkoIIYQtSUAJIYSwJZfVBUxFdna2XrJkidVlCCGECIH6+vozWuuc4ffPyIBasmQJdXV1VpchhBAiBJRSh0a6X7r4hBBC2JIElBBCCFuSgBJCCGFLElBCCCFsSQJKCCGELUlACSGEsKUZOc1cCDG2rq4ufD4fra2t9PX1WV2OiEFOp5O0tDSysrJISEiY0jEkoISIMl1dXRw+fJjMzEyWLFlCXFwcSimryxIxRGtNT08PLS0tHD58mEWLFk0ppKSLLwxONHdQ/tM6Trd2Wl2KiEE+n4/MzEyys7OJj4+XcBIRp5QiPj6e7OxsMjMz8fl8UzqOBFQYfO3RN/nLnlO86J3aP4oQ09Ha2kp6errVZQgBQHp6Oq2trVP6XgmoEHv7dCu/aTgKwJGmcxZXI2JRX18fcXFxVpchBABxcXFTHgeVMagQ+9qjb5Icb97WI74Oi6sRsUq69YRdTOd3UVpQIbT7qJ9HXjvJP1y/FE9OCkelBSWEEFMmARVCX/nLG2Qmx/HhdUtZkJnE0SZpQQkhxFRJQIXI8/vP8MxbZ/inm5eRlhjHwsxkjjV10N+vrS5NiBmtoKCAqqqqUR/Py8ujoqIiZMcbS0VFBVVVVVRXV0/6e6fzurFKAipEvv7om8zLSOSetYsBWJCZRHdfP6dbuyyuTIjotnXrVsrLy8P+OrW1tXi9XsrKysL+WsKQSRIh8OrRZuoONfH5d19KYpwTgAVZyYCZyTc3I9HK8oSIakVFRRF5Ha/XS1ZWFgButzsirxnrJKBC4MHnD5Ic76S4cMHAfQszTUAdbTrHlUuyrCpNiAFffHgPe4+3WFrDpbnp/MftKy2tYaqysrIGdvIuLCy0uJrYIF1803S2rYuHdx/nrvwFpCeeP/dkQWYSIFPNhQiF/fv3U1JSQmZmJgUFBXi93oHHho/teL1eCgoKUEqxYcMGSkpKyMvLY8uWLRM63miKiorwer34/f4hLajy8nKUUpSXl0/oOEF+v5/y8nIyMzOHjKOVlJQM6bLcsmULmZmZQ743MzOThoaGUY9dVVVFeXk5fr+f6upqqqurKS8vH/N77EhaUNP0y51H6O7t5++uXTzk/sQ4JzlpCTLVXNjGTG25gPnAPXDgAG63m4KCAioqKtixY8eIzy0pKWHjxo3U19dTUVFBQ0MD+/fvn/LxBissLGTz5s1UVlYO3FdZWUlBQcGkx6bWr19PUVERTU1NA3WXl5ezcePGIZM+ampqKCwspKGhgfz8fGprawHIz88f8bi1tbWUlpZSX19PSUkJNTU1A49t3rx5Qj+nXdimBaWUyldKFSuliq2uZaJ6+/p56IVDrFuWzbLZaRc8viAzSVpQQoRAaWnpQKtl48aN+P3+UZ/b0NAwEBYbN24c6Jab6vHAtHa8Xi87duxgy5YtQ1pKtbW1kw6n2tpa/H7/kKDbtm0bVVVVAy214Gt4vV5KSkrYvn07YAKrtLR01GN7PB7cbjd1dXVDjj94DG2msE1AAZ/VWlcDWUopj9XFTETN3lOcaO7k765dMuLjCzOTZbkjIUIgLy9vyO2xFh8tKioa6PLbvn37iJMoJnM8MK2b/Px83G43mzZtGuiCa2homNIkjYaGBjyeoR9zwcCsq6ujqKiI6upqamtrKSoqorCwcGBqe3V1NSUlJaMeO3hcr9c7pJVVU1PDhg0bJl2rlWwRUEqpMmCnUsqjta7SWk+8I9dCD714iAWZSdyyfPaIjy/MSuJEcye9ff0RrkyI2FZTU0NeXh5er5dt27ZN61i1tbUD3WpguvR8Ph8VFRX4fL4pzejzeDwXtOyCrbjCwsKBrrlgqOTn5+Pz+WhoaMDr9Y4birW1tRdM5AiG3Uxii4AC8oBZgE8ptVUp5R7+BKVUmVKqTilV19jYGPEChzvd0snz+89yV/4CnI6R15pakJlMX7/mRLNsuyFEpHi9XrZu3Up9fT07duyY9pTwYItk8ASDyspKtmzZMqRrMNhtNxHFxcVkZWUNtMT8fj8lJSUUFxfjdrspLS2lrq6O2tpaiovNqEdpaSmbN28euD2W4a2l4PjVTJseb5eAAtivtfYD9cAFHbqBllWh1rowJycn4sUN98fdJ9Aabl+TO+pzzk81l3EoISLF7XaTl5dHZmYmSikyMzMntdLEcB6Phx07dlBRUcGWLVuoqqrC7/fT1NRETU0N5eXlVFRUXDC7bzz19fUDe3ctXbqU/Pz8gQkMbrcbj8czZMxow4YNVFdXs3HjxnGP7fV6hwTZaF2dtqe1tvwLKAbKAtc3Ba+P9lVQUKCt9t7vPKtv++bTYz7nQGObXlzxR7195+EIVSWE1nv37rW6BMts3bpV5+fn66ampoH76uvrNaDr6+utK8xi+fn5uqamxrLXH+93EqjTI3zW26IFpc3kCLdSqihw29YLVh3xnePlw/4xW08Aue4klJIWlBCRNnzSg9vtnnHdW6E01ckcVrPNeVBa6+BZdLVjPtEGHt59HIB3XzZvzOfFuxzMS0/kqE9m8gkRCcHp3iUlJQPTtAsLC9mxY8cFs+ZiRUNDw4TGrezINgE1kzz8ygnyF7lZGFhvbywLZKq5EBFVVlYmC7oOMnhsa6axRRffTPL26Tb2nWgZt3svaEGW7AslhBBTIQE1SY/tOwXAravmTuj5CzKTOdnSSVdvXzjLEkKIqCMBNUlPvdnI8rlpzMtImtDzF2YmoTWc8Mu5UEIIMRkSUJPQ3tVL3cEmbrh44udhLRy0L5QQQoiJk4CahBe8Z+nu6+eGiyYeULLthhBCTI0E1CQ8/WYjSXFOCpdkjv/kgHkZSbgcSrbdEEKISZKAmoSn3mxkrSdrYFv3iXA6FLnuJI7ITD4hhJgUCagJOnS2nYNnz3HjJMafghZkJkkLSohpqq6uHlhjL7j1xFiG77Q7/PZwg3e1DVdNY6moqKCqqmpKxxnvZ5up5ETdCXr6rTMAk5ogEbQwM5nHXj8d6pKEiBler5f77rtvYBfccNi6deukVpsIZU21tbV4vV4qKyunHXTRRAJqgnYe8DEnPYGl2SmT/t4FmUmcaeuio7uPpPiJdw8KIYzg/kbhXE9vsmvVhbKmwbvdxvKagcNJF98E1R9qonBxFkqNvPfTWBbNMlPND8uafEKIEWRlZQ1sYDh8o8FYJi2oCTjV0skxfwcfvG7JlL7/4jlpALx+soVL5qaFsDIhJuGRz8DJV62tYe5quO3Lk/qW4D5MAJmZmVRWVlJWVkZBQQEbN25k06ZNAGzZsoXt27dTX18/pdIKCgooLy8fWMcvePyamhrq6uoG9oXyeDyj1uT3+6moqOBXv/oVWVlZFBcXU1lZOe5rFxUVcd99912wp1R5eTlVVVWUlZVRUVEx4S7I0eooKSkhKyuLrVu3AuY927x5M01NTQPfm5mZyWOPPTZku/jBqqqqqK+vp7KycmCn4eC+WKN9z1RJC2oCGg6Zf7yCxROfXj7YstmpxDsd7D3eEsqyhIgJlZWVbN26lfz8fJqamiK6EOzmzZvZsWPHwAd4cBLFaDWtX78et9tNU1MT+/fvx+v1DuyaO57CwkI2b9485L7g60x2fGy0OjZu3Dhk+/qamhoKCwsHdgsOPjZa0NTW1lJaWgowsANwcXExGzZsuKD2UJAW1ATUH2oi3uVgZW7GlL4/zung4rmp7JGAElaaZMtFmG3Wgy2aYGtqNMEt3we3mLZt2zbQwhptbMnv9+P1etmxYweZmZmUl5cPhFFtbe2kA3m8OrxeL16vF4/Hg9frpaKigu3bt5Ofn09NTc1AAI3E4/Hgdrupq6tj27ZtA/cPHkMLJWlBTUD94SYum59BvGvqb9fKeRnsOd4c3EFYCDED5OXlDbk9fCPEwRoaGi5o5QRDKTi+NJKSkhLy8/Nxu91s2rRpoMU11U0Gx6ujqKiI6upqamtrKSoqorCwcGDmYHV1NSUlJaMeO3hcr9c7pJVVU1PDhg0bJl3reCSgxtHZ08drx5qn3L0XdGluOk3nejjZIovGChGNPB7PBUHk9/uB0Sc+1NbWDulyq6ysxOfzUVFRgc/nm9KMvvHqKCkpoaamZiBU8vPz8fl8NDQ04PV6xw3F4OzF4feFY8deCahx7DneTE+fJn+aAbUyN90c75h08wkRCh6Ph7NnzwLmA3j79u2W1lNcXExWVtZAC8jv9w+M04wWNMEWSXAMCExIbdmyZSBU4Hy3XSjqKC0tpa6ujtra2oGddktLS9m8efOEdt4d3lpqaGgYaAGGmgTUOHYdaQbgioXuaR1n+bx0lIK9JySghAiF4Ay3vLw8SkpKbDE9u76+Hp/PR2ZmJkuXLh13N9vgzMDgrMCqqir8fj9NTU0DM+MqKioumN03nTrcbjcej2fImNGGDRuorq5m48aN4x7b6/UOCbLt27eHpfUEoGbimEhhYaEeq083lP5txys88cZp6u6ffv/qzV99kovnpLL1Xuv/I4notW/fPlasWGF1GSJGFBQUUFlZOWZIjfc7qZSq11pf8MEoLahx7DvZwvK56SE51qW56dKCEkJElalO5pgICagx9Pb18+apNlbMC83JtZfOS+eIr4Pmjp6QHE8IIazU0NAwoXGrqbJdQCmlxj/tOkIOnm2nu7c/ZC2o4EQJOWFXCBENxhtjmy5bBZRSqgiY+OnSYbb3RCsAK+aFrovPHFcCSgghxmObgFJKeQDvGI+XKaXqlFJ1jY2NEanp9RMtuByKvNmTX8F8JLPTEslJS2DP8eaQHE8IIaKZbQIK8GitRw0orXWV1rpQa12YkzP5PZmm4vWTreTlpJLgCt0WGZfOS5cuPhF2M3F2rohO0/ldtEVAKaWKtNa14z8zsvadaAnZBImglbnpvH26ja7evpAeV4ig+Ph4Ojo6rC5DCAA6OjpISEiY0vfaIqAAn1KqSClVDHiUUqFds30KWjt7ONHcycUh3h5jZW4Gvf2at061hfS4QgRlZ2dz9OhRfD4fPT090poSEae1pqenB5/Px9GjR5k1a9aUjmOL1cy11g1gxpkAt7XVGAfPmM0FPdmpIT1ucKLEnuPNrJo/tdXRhRhLRkYGCQkJNDY2cvbsWXp7e60uScQgl8tFYmIiixYtIjExcWrHCHFN06K1rgKqrK4DwHvGtHCmssX7WBZnJZMS75RxKBFWiYmJLFy40OoyhJgWu3Tx2c7BM+dQChYHtmsPFYdDsWJeuuwNJYQQ45CAGsXBs+3kZiSRGBe6GXxBK3PT2Xeihf5+GRsQQojRSECNwnumnSXZoW09BV2am057dx+HfOfCcnwhhIgGElCjOHimPeTjT0HBrePlhF0hhBidBNQImtq7ae7oYcms8ATURXNScTmUjEMJIcQYJKBGcOBsOxD6GXxBCS4nK+dn8NIBX1iOL4QQ0UACagQHGsMbUAA3XpzDy4eb8J/rDttrCCHETCYBNYKjTWaZmPmZSWF7jZsuyaFfw9NvnQnbawghxEwmATWC4/4OZqclhHSR2OHWLHCTmRzHk2+cDttrCCHETCYBNYLjzR3Mc4ev9QTgdChuuDiHp95olPOhhBBiBBJQIzjm72C+e2prR03GTZfkcLa9m9dkurkQQlxAAmoYrTUn/J3kZoS3BQVww0U5KAVPvB6ZDRiFEGImkYAaxn+uh46ePnLD3MUHMCs1gcsWuHnyTRmHEkKI4SSghjnmNzP4IhFQADddnMOuI3587TLdXAghBpOAGuZ4IKDmRyigbl4+G63hmbekm08IIQaTgBrm+EALKvyTJAAum59BVko8T74hASWEEINJQA1zvLmTBJeDrJT4iLyew6G44aJsnnpTppsLIcRgElDDnGjuZG5GIkqpiL3mzctn42vvZvcxmW4uhBBBElDDNLZ2MjstIaKveX1gurmsKiGEEOdJQA3T2NpFToQDKislnjUL3Dwh41BCCDFAAmqYxtYuclIjG1AAN18ym91H/Zxt64r4awshhB3ZIqCUUm6lVL5SqlgpVWlVHZ09fbR09ka8BQVm2SMz3VxWNxdCCLBJQAGlQKHWuhpAKVVmRRFnAq0XKwJq9fwMZqXEU7PvVMRfWwgh7MgWAaW1rtJaVwVuegDv8OcopcqUUnVKqbrGxvCM1TS2WhdQDofi9jW5PLrnJKdbOiP++kIIYTe2CKggpZQH8Gmta4c/FgixQq11YU5OTlhefyCgUiNzku5wf3/tEnr7NQ+9cMiS1xdCCDuxVUABxVrrcqtevNHCLj6AJdkprF8+h4dePExnT58lNQghhF3YJqCUUsVa6y2B6/lW1BBsQc1KjcwqEiP50Lol+Nq7+cOu45bVIIQQdmCLgFJKFQGVSql6pVQ9kGVFHY2tXWSlxBPntO5tucYzi+Vz0/jRcwfQWpY+EkLELlsElNa6Vmudp7UuCHxdMAYVCWfaupgVoTX4RqOU4kPXLeX1k638df9ZS2sRQggr2SKg7KLpXA+ZFgcUwB2X55KVEs+PnjtgdSlCCGEZCahBmtq7yUyOs7oMEuOc3HP1Ih57/TQHz7RbXY4QQlhCAmqQpnM9EdtmYzz3rF2My6F48PmDVpcihBCWkIAK0FrjP9eNO9keATU7PZHbL8tlR90RWjp7rC5HCCEiTgIqoLWrl95+bYsuvqAPXreU9u4+frXziNWlCCFExElABfjbTSsl0yYtKIDVCzK4ckkmDz5/kD7ZbVcIEWMkoAKaznUD9googA9dt5SjTR088toJq0sRQoiIkoAK8AUDKsU+XXwAGy6dw/K5aXzx4b34AzUKIUQskIAKCH7422WSRJDL6eBrpWtoau/mP/6wx+pyhBAiYiSgApoCY1BZoQiocz54tRr6Q7Pg68rcDP75lov4/a7j/Fm6+oQQMUICKqDpXDdKQXrSNLr4utvh6a/At9bArz8M+58IWX0fvTmPVfPT+fffvibbwgshYoIEVEBzRw9pCS6cDjX5b+7rgZe2wbcuh8f/D8wPLMbuD92+TnFOB18ruZzWzl7+9+9fk4VkhRBRTwIqoLWzd/Ktp/5+05X3wJXwp3+F7IvgwzVwz29AOaHlWEhrvGRuGp/YcBF/evUkD++Wrj4hRHRzWV2AXbR09JCeOMGA0hr2Pwa1X4STu2HOari7GpYVgQq0wNLmQXNoAwqg7HoPf9lzis///jXWerKYnWbN7r9CCBFu0oIKaOnsIT1pAnl9tA5+fDs8dBd0NsOd26D8abhow/lwAsiYH/IWFARm9ZWsoaO7j8/9Rrr6hBDRSwIqoLWzl7SxWlCNb8Av74YfrIfG1+G2r8DH6uCyUnCM8DamhyegAJbNTuXf3nEJtftO8ZuG8LyGEEJYTbr4Akbt4ms+Ck9+GXb9DOJS4OZ/h7UfhYTUsQ+YMR/e+JPpDlRTmHgxjg9et5S/7DnJ/b97jbkZiVy3LDvkryGEEFaacAtKKXWLUuryYfctVUrdqZRaEurCIq2ls3doF985Hzx6P3w7H3Zvh6v/ET6+C27cNH44AaQvgN5OOBeeXXGdDsV37s5nUVYyH3xwJ7V7T4XldYQQwirjBlQgmPqBGqBBKdWnlPo0gNb6APA4sD+8ZYZXX7+mrSvQxdfdDk9/1UwZf/4BWHUX/HM93LoZUibRSsmYby6bj4alZoDZaYlsL1/LirlplD9Uz+93SXefECJ6jBlQSqmlQBVQDiwDCoCNwNVKqbeUUjdrrf1A6PuwIqitsxcXvVx95jfw7Svg8f+CJdfBR56H930P3Ismf9D0QECFaRwqyJ0cz8/uW0vh4kw+sX0XP3/xcFhfTwghImW8MagyoEBr3TzovpeBagCl1F1KqfuA2jDVF379/fTs3kFt/BdY8sYpWHQNlP4EFq2d3nGDARWGqebDpSa4+PGHruIfH6rnc799lfauXu67wRP21xVCiHAaL6DqhoXTEFrrXweubptuIUqpYsAPeLTWVdM93ri0hv2PQ+0XyD65m0YWUn/d9yko+pvQTGpIyQFHHLSEr4tvsMQ4J1X3FvLJ7bv40p/20dbVyyeKLkKFYYKGEEJEghrrPBql1J1a69+EvQgTTmitq5VSZYBXaz1qq6wwLU3XFRRM+fXq56WzYP4bzEk9SVt3Ko/7r+MT6ffy0N4dXNsSwi6yq/dASwrsWxK6Y46jD8VnPO9gx+zVfOhEHZ879CQu5FwpIYTRi6LL4aLb4aRLmcte5aBPOegJXPYO+urDQY/DXA7cpxz0KViY5WXWSReelvZp1aSeeqpea104/P7xWlBZ03rVibsS2B647gXyGdZtGAiuMoDLEhKm9WJJSw7jjD/HF3o+wM/719OdHodD97Oga9TG4tR0xUNCZPdwcqKp9P6Z1L4ufjSvkGczFvOFg49xbYtsGy+EXWigRznocMTR6XDR4YijwxlHh8NFpyOOc444OpyuwOPm/g5HHF0OpwkX5TQBM/h6IGy6hly/8PF+Nb3TX5308T7ns/yL8zcscjRS7yyEltC8L8ON14L6MvALrfUrIzx2C7ABEyjf11o7p1yEUluBrVrrBqVUEbBBa10x2vMLCwt1XV3dVF+OpmNv0RnnRsefny6eHO8M/V5Qv74PjrwAn3g1tMedAK01j+49xX/9cS9Hmzp41+p5fO5dK5jvTop4LULMRFprunr7aevqpb2rl/auPtq7L7ze0dNHZ3cfHT19nAtcdvb00RG43tHTP/B4x6D7+/on37MR51QkuJzEuxzEOx0kxA27DD7mcpAwcOkkYdDtod/nJM7pIM6pcDkcOB0Kl0PhGnzbGbgPTdaBP5DT8A3imw/SNfsyWq/ZRNKKW0mZ6DJxo1BKTakFtRk4oJR6FNOicQN5QGng9n1a65ZAwEyHn8i11sicf1FkXihjPuw5YRaVHWm1iTBSSvGOlXO58eIcqp728t0n3+ax10/x0ZuWUXaDh8S4Kf89IYQt9ffrQGgMC5Ku3vP3D1zvpb27b4TnnP/+c90TDxGlIDnOSVK8k8Q4J0mDrruT4khKTxz2mIPkeNeQ20lxQ7934Pag61PabWG6+vth72/NggVn3oQ5q+C2n5NwyTtJCPMY95gBpbVuVkoVAjswoQTmfKgirfXLYGbyYQJmOnZiwg/AE3iNmS99PvT3QPtpSJtrSQmJcU7+Zf1F3Jk/n81/ep2v17zJjvoj3P+uS/lfl86RSRTCFvr6NW2dvbR09tDa2Utr8LIreHvwY4MeH3RfW1fvhF8vKc5JSoKLlAQnKfEuUhNcZKXEszArmZT4wGPxLlISXKQmOEkOXE9JGPyY+d6keNNCibr/S1rDvofhyc1wei/kLIeSH8OKOyL2B/e4Sx1prb2Y858uoJTKAGq11tNq/QQmR2wKdO+5x5ogMaMMnmpuUUAFLchM5jt353P322f4wsN7KP9pPeuWZfOBaxZz4yU5JLikRSWmprevf4QQGR405nrLsIBpC1xv7x5/9+l4l4P0RBdpiXGkJbpIS3SRk5pKWqKL1OD9Ca4hwTNSqCTHT3Hft1ihNbz5Z3jiv81uDbOWwV0/hJXvA0dkPyemtRZfcAq6Uipdaz2tYTKt9ZbA1egIJzi/mkTLMUbJ+Ii7dlk2/+9fruehFw7xwONvU/bTetITXbxz9TzuuDyXq5fOkv+8MUZrTWdPP80dPTR39NDS2UPzuZ6B24Pvbxl+X4cZgxlPgstBWmJcIGBMmMxNTxy4PvgyfYT70hJd8kdUuGkNbz8GT3wJjjdA5hJ47/dhdQk4rVm2dVqvqpRKB7YA/zDdY0Wl9AXmMsyrSUxWnNPBB69byj1rF/Pc22f4/a7j/OGV4/xy5xHmpidy+5p5vOfy+azMTY++bosopbVZrquls3dIuLQMDp2OEUKno5eWjh66+/rHPH5agov0pDgyAl9Ls1PISIojPTGO9KQLwyR9UMCkJriId8nGCbalNRx4yrSYjrwIGYvgjv+BNe8H5/QmP0zXlEIlEEyfBTYBB5jhSx2FTXIWuBLDuh7fdMQ5Hdx0yWxuumQ2Hd191O47xe93HefB5w+y7ZkDeHJSeNfqeVyxyM2q3Axmp8vmiOEU7CoLhklLR++o4dLSMbQ109LZO+aAvkMxEDDpieYyNyNpSOikJ7kGrg/+Sk1w4XJKwESlg8+ZFtOh58yQxLu/AZffA64Qz2ieokkF1AjBVAo8BoRnye6ZTqmw7gsVSknxTm5fk8vta3Lxn+vmkddO8vtdx3jgibcJnokwOy2BVfMzWJWbbi7nZzAvI1FaWQFaa9q7+863WgaFR8uQbrLeIV1mLYHnjDfI73KoQWEShzs5nsWzUkYMlyFhlBxHarwLh3TdiqAjL8Hj/8e0nFLnmv3t8j8Acfb6I3RCATUomCowK5eXBpc5CkyUEKPJmB+R9fhCyZ0cz/uvWsT7r1pEW1cve4+38NqxZl473syeYy08+cZpgn+sZ6XEszI3nWWzU5mdlkhOWgKz0xLICXxlJcfPiA/G3r5+2rv6aO3qoa3LDN63dfWOen3EoBmnFQPnu8rSk8x4zKKs5MD1862Y89eHBk9SnFP+GBDTc7QenvxveLvWLMf2jv+Gwg9BnD3PjxwzoALB9DlMi6kecwLtY5EoLGqkLzB/pcxQqQkurlqaxVVLz0/U7OjuY9/JFvYca+a1Yy28eqyZlw8fHbEF4HQoslPjTWClJpCdmkByvHPgBMLBJxSOdJ8C+rWmrx/6+vvNpdb092t6+81ln9b09Zuv7t7+8ydKDjphsrOnf+BEyYH7es1lW1cvnT1jj8EEJQemIGcEAiY7NR5PTsqQgBk8NnP+uhmPkQkowhInXoEnNsObj0BSFhR9Ea66D+JTrK5sTOO1oB7HrMohwTRV6bnQehL6ei2bCRNqSfFO8hdlkr8oc8j957p7aWztGvg6Peh6Y1sXp1s72Xeilc7ePrp7++nq7Z/S2fQTEedUAyc9DlzGO0mKczArJZ6kTCeJLnNfcGpyaoKZrpw66PrAY4lmmrIEjJhRTu0x5zHtexgSM+CW+83mqwlpVlc2IeN9Yq4HikBWG52yjPmg+6Dt1Plp51EqOd7F4lkuFs+a+F9lvX39dPf1DwSWueyjK3BbawaWX3EohdOhcDrA6XDgVAqHg8B9CqdSxLscJAaWbxEiZjW+YYJpz28hIR1u/Ays/Qgkua2ubFLGXUkC+LVSKiOwYkST1vrxyJQWJQZPNY/ygJoKl9OBy+kg1MsgChGTzu6Hpyrh1R3gSoLrPw3XfMzMKJ6BJtTnNCyo7gP2S1BN0OCt3xdeZW0tQojodHY/PP0V2L0dnAkmlK77OKRkW13ZtExqUCQQVNsGBdVZzDRzMZoIbf0uhIhBPi88/VV45ZfmpNqrP2KCKW2O1ZWFxJRG7QcHFWaPJhk5Hk1iBsSnzrip5kIIG2s6aFpMu34BDhdcVQbrPmH5mp+hFoq1+L6ilAr/Fu0z1cDJuvZcTUIIMYM0HYJnvgq7fg7KCVf+A6z7JKTPs7qysAjJvOfgorFiFBnzoeW41VUIIWYq/2HTlbfrZ6Ac5uTadZ80p7FEseg4Mcfu0nPh1F6rqxBCzDT+I/DM1+Dlh0xvTMHfw7pPxcyMYAmoSEhfYM6D6u22zSKMQggbaz5mgqnhJ+Z2/r1mynjGAmvrijAJqEjImA9oaD0BmYutrkYIYVctx+GZr0PDj0H3wxX3mGByL7K6MktIQEVClsdcnn1bAkoIcaGWE/DsN6D+QbPyzOV3m2CK8c8LCahIyFluLhtfh2Xrra1FCGEfrSdNMNX9X+jvhcv/Fm74V7ObrZCAioiUbEjONgElhBCtp+C5b0Ldj6Cvx+xee8O/QtZSqyuzFQmoSMlZDqcloISIaW2n4blvwc4fQl8XXPY3Jphm5VldmS1JQEXK7OWwewdobaaLCiFiR1sjPP8teOkHgWDaCDf8mwTTOGwRUEopN+AJfF2pta6wtqIwyFkOXc1mJl+Un1wnhAhoPxNoMf0AejthdQncsAmyl1ld2Yxgi4ACSgG01lVKqSuVUmVa6+haPmnwRAkJKCGiW/tZeP7b8NI26DkHq4tNMOVcbHVlM4otAmpYGHmAmuHPUUqVYRamZdGiGXhOwOwV5vL065B3i7W1CCHC45wPnv8feKkKutth1Z1wYwXkXGJ1ZTOSLQIqSCnlAXxa69rhjwVCrAqgsLBw5u3wm5INybOgcZ/VlQghQu2cD/76HXhxK3S3wcr3wY2bzv9hKqYkYgGllCoGhm/r6B0WRsVa6/JI1RRxOSvMVsxCiOjQ0QR//S68+H3oaoFL32taTHMutbqyqBCxgNJaV4/1uFKqWGu9JXA9X2vdEJnKIijnEni1WmbyCTHTnfPBC987H0wr7oCbPgNzVlpdWVSxRRefUqoIqFRKfTZwV/TN4gPT3O9qNmePR+n+LUJEtSFdea2w4nbTYpq72urKopItAirQzRf9JwQEB0ob90lACTGTtJ+Fvz4QmPzQBpe+x8zKm7vK6sqimi0CKmbkBAZMG9+QmXxCzATtZwKz8gLTxVe+z5xgK2NMESEBFUmpOWYm32mZySeErQVXftj5Q+jpgFV3mWCavdzqymKKBFSk5SyXRWOFsKvWU+YE2+BaeauKTTDJCbaWkICKtJzl8JrM5BPCVlpPmiWJ6n4Efd2wutQs4pp9kdWVxTQJqEjLWQ6dMpNPCFtoOWG2vah/0Gx7cdlGWV3cRiSgIm324DX5JKCEsETzsUAw/TiwUeD7zQ62wd2vhS1IQEXa4EVj8262thYhYk3zUbODbcNPQPebHWyv/7TsYGtTElCRlpIDSVkyUUKISPIfgWe/Dg0/NbevuBvWfQoyF1tblxiTBFSkKSW76woRKU2HTDC9/DNzO/9eWPdJcM/AHRFikASUFWYvh9d+LTP5hAgX3wF45mvwyi9AOaDg70wwZSywujIxCRJQVshZYWbytZ2CtLlWVyNE9PB54elAMDlcUPghuO4TkDHf6srEFEhAWSG4Jt/pfRJQQoTCmbdNi2n3dnDGwVVlcN3HZabsDCcBZYXZg9fkk5l8QkzZqb3wzFdhz2/BmQBXl5tgkj/8ooIElBVSciApE07vsboSIWamE6/A01+BfQ9DfCpc+y9wzcfMepciakhAWUEpWHAlHH7B6kqEmFmO1sFTW+Ctv0BChtnyYu1HIHn4Zt0iGkhAWWXJOnjrUbM4Zdocq6sRwt4OPmdaTN4nTO/DLfebcabEDKsrE2EkAWWVJevM5aHnYNWd1tYihB1pDd4nTTAdes50jW/4Tyj8MCSkWl2diAAJKKvMXQPxaXDwWQkoIQbT2vQuPP0VOLoT0nLh1krI/wDEJ1tdnYggCSirOF2waK0JKCEE9PfDG//PBNOJVyBjEbzr63DFPeBKsLo6YQEJKCstWQe1/wFtpyF1ttXVCGGN/j4zTfyZr8HpvWZF8fd8x2x94YyzujphIQkoKy253lweeg5Wvs/aWoSItL5eeHWHCaazb0H2JXDnNlh5p+lhEDHPdr8FSqlKrXWF1XVExLw15hyOg89KQInY0dsNr/wcnvk6+A/BnNVQ8mNYcQc4HFZXJ2zEVgGllCoCYmfHMBmHErGkpwNefgie/Sa0HIXcK+DWL8Mlt8miyWJEtgkopZQH8I7xeBlQBrBoURQtlb9kHdR+Adoa5Sx4EZ06m2HnD+GF70J7IyxcC3d8C/LWSzCJMdkmoACP1rpWjfILq7WuAqoACgsLdSQLC6uBcSjp5hNRpv0MvPA9eGkbdDWbQLr+07D4WgkmMSERCyilVDEwfD0SbyCUirTWtZGqxVbmrYG4FHOmvASUiAbNR+H5/4H6H0NvJ1x6h9mLKfcKqysTM0zEAkprXT3Gw77A+JMb8Cil8rXWDZGpzGLOOBmHEtHhzFtmfGn3L83tyzaavZhyLrayKjGD2aKLLxhGgXEmt7XVWGDJOnjsizIOJWamE6+YGXl7f29OqC38MFz7MdlWXUybLQIqaPA4U0wZcj7Uey0tRYgJO/S8OYfp7VpISDfdeGs/Kn9kiZCxVUDFrNzLzTiUBJSwO61NID3zNTj8V0jOhvWfhyv/QVYWFyEnAWUHzjhYdLWMQwn76u8zXXjPfh1OvgrpC+C2LXDFvbKAqwgbCSi7WLIOHvtPMzU3JdvqaoQwervNpIdnvwm+/TDrInjPd2F1Cbjira5ORDkJKLsYPA516XusrUWI7nZo+ImZLt5yzJwOUfJjWHE7OJxWVydihASUXeReAXHJpptPAkpYpaMJXvoBvPg9OHcWFl8Hd3xbVn0QlpCAsgtnHCy9EfY9bNYnk79SRSS1njJLEe38IXS3wkXvgOs/Zc7RE8IiElB2cvnfwpuPwP7H4aINVlcjYsGZt+H5b8Mrv4D+Xrj0vWa6+LzLrK5MCAkoW7n4VkieZVZ8loAS4XS0Dp77Juz7Izjjza6113wMZuVZXZkQAySg7MQVD6tLoe6HcM4HycOXLhRiGrSGt2rguW+ZxYkTM8zirVeXy47OwpYkoOzmirvNAPWr1XB1mdXViGjQ1wOv/doE0+m95hymd2yG/HshIc3q6oQYlQSU3cxdDXMvg5d/KgElpqerzUwV/+t3zAaBsy+F922FVXeZSTlC2JwElB1dcQ88ssmcsT93tdXViJmmrRFe/D7s/AF0+mHxOnj3N8y4pkwVFzOIBJQdrS6BR++Hl38Gt33Z6mrETNH4Bvz1AXhlO/R1w4p3m+0uFhRaXZkQUyIBZUfJWXDJO+HVX8GG/5QlZcTotIYDT5tgeutRcCWaccy1/wTZy6yuTohpkYCyqyvugb2/gzf/bHYkFWKwvh7Y81uzFNHJ3WZV8Zs+B1d+WNZyFFFDAsqu8m6BtHnmnCgJKBHU2Qz1D8KLW80aedmXwO3fNrvXxiVaXZ0QISUBZVcOJ6z5GzM1uPUkpM21uiJhpaZDZuJDw0+gu80sLvzub8CyDeBwWF2dEGEhAWVnl98Dz34DXvklrPuE1dUIKxyrh+cfMHsxgZkifs0/mU0uhYhyElB2lr0MFq6FXT+D6z4uU4RjRV8v7PuDaTEdedFsp37NR+Hqf4SMBVZXJ0TESEDZ3RV3wx/+GY7uhIVXWV2NCKdzPjO+tPMHZnwpc4lZ8eGKeyAx3erqhIg4CSi7W/k++Mv98MSX4N7fSSsqGp3eBy98D3b/Cno7YOkN8M6vwsXvkG1XREyzTUAppfIBD4DWutricuwjIQ1u+XezssS+h2VGX7To7zfnLb3wXTjwlDl/6bJS0403Z6XV1QlhC7YJKOCzWusSpVSZUsqjtfZaXZBtFH4Y6n8Mf/l3WFYE8clWVySmqrMFdv0cXtoKPi+k5cL6z0P+30PKLKurE8JWbBFQSqkyYGcgmKqsrsd2nC545xZ48F1mD5+bP2d1RWKyzu6Hl7aZ89q6W2HBVXDL/bDiDlm4VYhR2CKggOAuaT6l1FagQmvtH/yEQIiVASxatCiy1dnBknWwqhie/SaseT9kLbW6IjGevl546y9mG/X9j4HDBSvvhLX/CPMLrK5OCNtTWuvIvJBSxcDwHfi8WutapVQlsF9rXRUIIrfWestoxyosLNR1dXXhLNeeWo7D/xSC50Z4/y+srkaMpvWUOaG2/kGzzUXaPCj4IOR/ANLnWV2dELajlKrXWl+wqnHEWlDjTHzYyfnwcgP+cNczI6Xnwo3/BrVfMDujyrbw9qE1HHrOTBHf9zD094LnJrMa/cW3mW5aIcSk2OJ/jda6Wim1SSlVFLgt41CjWftRM47xSIWZjuxKsLqi2NbZbLa3qPshNL4OiW4zE6/gg7KauBDTZIuAAhjUpVdraSF250qA2yrhobvMTqnXf8rqimLTid2mtfRqNfS0Q24+vOe7sOpOiEuyujohooJtAkpMwrIiuORd8PRXzCrWGfOtrig2dLWaLS4afmJW9nAlweq7zGkA8/Otrk6IqCMBNVPd+t/wwFXw589A6U9khYlw0RoO/9V0q+75LfScg+yLzRJEl78fkjKtrlCIqCUBNVNlLoGbPgOPfdGcwPuOL0lIhVLLcXjlFyaYfF6IT4PVJXDFvWYLdXmvhQg7CaiZbN0noe0UvPAdMza1/vPywTkdvV3wxiNm9fi3a0H3w+J1cMMms8RUfIrVFQoRUySgZjKl4NYvQ28nPPt1Mzh/4yarq5p5Tr5mWkq7t0OHD9Lnw7pPweV/C7Pyxv9+IURYSEDNdErBu75h/vp/4kumJXXdx62uyv6aj5kxpVd3wIld4IyH5e8yW1t4bpZVxIWwAQmoaOBwwB0PmJCq+bxZGfvqcqursp+202Zn2td+A4efN/fNWwO3bTHjS8nDFzoRQlhJAipaOF1wZxX0dZutOZzxUPhBq6uy3jmfWdnhtV/DwWfMuFLOCrj5fnPOknThCWFbElDRxBkHxT+CX94Nf/ykaUld/n6rq4q8zhZ4408mlPY/bpYdyvLA9Z82i7XOudTqCoUQEyABFW1cCbDxp/DzjfC7j5iurFs+D6k5VlcWXud88PZjsPd3Zp3Cvi7IWGiWhlp1l+nKkxmOQswoElDRKC4J3v9LM2nixe/Dnt/DTRVwVVn07D3U3w8nX4G3as3OtMfqTPdd6hzTtbnqLphfaMbnhBAzUsS22wilmN1uYyoa3zSrTex/zKyAcOuXYdl6q6uamo4m2P+EaSG9XQvtp839uVfARf8Llm0wSw7JDDwhZhTLt9sQFsm5GO75Nbz5Z/jzZ+GhO+GSd5qVJ7I8Vlc3Nq3h5KumhfR2LRx5CXSfWTF82XoTSMvWQ+psqysVQoSBtKBiSW8XvPBdePqrZrbf2o+Y6dVzVtljfKb9DBxrgOMNcKzefJ07ax6bt2ZQK6lA9lcSIoqM1oKSgIpFLSfMGn6vBHblTZltNtfLu9mcpBqJXV+72uDEK0PDyH848KCCnOUmiBZfa1ZvT5sT/pqEEJaQgBIXajluxnS8T4D3SWhvNPfnLDdBlXczLLoGEtIm38Lq7zPHaz0BrSfPXzYfheO7oHGfmdQAkLHIjB3NzzehNG+NeU0hREyQgBJj6++H03vMeUP7nzBbTPR2msccLohPNaERnwoJg68HLvt7hgZR26nzATRAmVl2c1eZIJpfYDb6i/Yp8EKIMckkCTE2hwPmrjZf130cejpNSB1/2WzU191muuW6W81lZ4tpgQXvc7ggbR6kzYU5K89fH3yZMlvGjoQQEyafFmJkcYmmiy/vZqsrEULEKDmLUQghhC1JQAkhhLAl23TxKaWKAT/g0VpXWVyOEEIIi9miBaWUKgK8WutawKuUyre6JiGEENayRUABdcCOQDB5tNYNVhckhBDCWrYIKK21H9gK7AAKRnqOUqpMKVWnlKprbGyMZHlCCCEsELETdQNjTMP31PZqrWsDXXw+rXWDUqoS2Km1rh7tWHKirhBCRA/LT9QdK3CAfK31lsD1zUBpBEoSQghhY3aZxVellCoDvMgsPiGEEMzQtfiUUo3AoWkeJhs4E4JyooW8H+fJezGUvB9DyfsxVCjej8Va6wsW5ZyRARUKSqm6kfo8Y5W8H+fJezGUvB9DyfsxVDjfD1vM4hNCCCGGk4ASQghhS7EcUDIRYyh5P86T92IoeT+GkvdjqLC9HzE7BiWEEMLeYrkFJYQQwsbsch5UWI23UnosraQ+1s+qlHIDnsDXlVrriogXGGET/bdXSlXK+wHB9TJh3JPvo4J8dpwX+FnLtdYbxnjcTwjfi6hvQQXeNAIrpQdXTp/w49FkAj9rKVAY/OAJnDwdtSb6bx+43xPB0iwxwffjs4HfjyylVFS/JxP47IipXRjG+oMkXJ+jUR9QwJWYFSoIXA7/JRrv8Wgy5s+qta4a9JePZ9Bzo9W4//aBD+Fofx+Cxnw/An+w7FRKeQK/K9H+voz3+yG7MJwXls/RWAgo97Dbsyb5eDRxD7s94s8a+FD2Bf8aimLuYbdHej88MfBBHOQednv4+5EXuM+nlNoa6BKOZu5ht4e8HxPZhSGGuIfdDsnnaCwElJ8LV1GfzOPRxM/EftZirXV5mGuxAz9jvB9KqaIYCOnB/Iz/+7E/8MFcD0R1FzAT+P0AarXWeYA/2M0Vo/yE4XM0FgJqJ+fT3QPUTPLxaDLuz6qUKg6uLB/tfeqM/374lFJFgQ8ej7wf7Bx03Y35UIpm470f+YO69TYTO3/ojiQsn6NRH1CBgT1P4K8d96BBvJqxHo9G470XgfsrlVL1Sql6ovw/3AR+NxoC92VxYRdG1Jng/xV3cAA82metjfd+ENiFIfB4abS/H4Gfs3BwSzHcn6Nyoq4QQghbivoWlBBCiJlJAkoIIYQtSUAJIYSwJQkoIYQQtiQBJYQQwpYkoIQQQtiSBJQQQghbkoASQghhSxJQQgghbEkCSgghhC1JQAlhI0qpYqXUfqVU07Cv/UqpSqvrEyKSYmLLdyFmguAOtVrrPKVUmda6KriCeoxvhidilLSghLAP36Bttd2ByyJiZ0dfIYaQgBLCJgIbAQ60pALygvcLEWskoISwnwoguJ+OZ6wnChHNJKCEsJ+iwWNOwQ0ChYg1ElBC2Ehgt9LqQXc1EAO7+QoxEtlRVwghhC1JC0oIIYQtSUAJIYSwJQkoIYQQtiQBJYQQwpYkoIQQQtiSBJQQQghbkoASQghhSxJQQgghbOn/A1sVouQ4EbOIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# To make figures of welfare comparison \n",
    "\n",
    "# linear interpolation \n",
    "index = next(i for i, x in enumerate(pi_array_3) if x > pi_array_0[1])\n",
    "print(index)\n",
    "\n",
    "Omega_array_3_new = np.zeros(len(pi_array_0)) #living on the grid of pi_array_0\n",
    "\n",
    "for j in range(1,len(pi_array_0)-1):\n",
    "     if pi_array_3[index-2+j] < pi_array_0[j] <=pi_array_3[index-1+j]:\n",
    "            Omega_array_3_new[j] = ( Omega_array_3[index-2+j]+(Omega_array_3[index-1+j] - Omega_array_3[index-2+j])\n",
    "                         *(pi_array_0[j]-pi_array_3[index-2+j])/(pi_array_3[index-1+j]-pi_array_3[index-2+j]))\n",
    "\n",
    "Omega_array_3_new[0] = np.nan             \n",
    "Omega_array_3_new[-1] = np.nan\n",
    "pi_array_0_new =list(pi_array_3[:index-1])+list(pi_array_0[1:])\n",
    "# Omega_array_3_new =list(Omega_array_3[:index-1])+list(Omega_array_3_new[1:])\n",
    "# Omega_array_0_new =list(np.zeros(index-1))+list(Omega_array_0[1:])\n",
    "\n",
    "diff = Omega_array_3_new - Omega_array_0\n",
    "Omega_array_00=np.zeros(index-1)\n",
    "Omega_array_00[:] = Omega_array_0[0]\n",
    "dif_00 = Omega_array_3[:index-1] - Omega_array_00\n",
    "diff_new = list(dif_00) + list(diff[1:])\n",
    "\n",
    "diff2 = aveOmega_0 - Omega_array_0\n",
    "\n",
    "plt.axhline(y=0, color='r', linestyle='-')\n",
    "plt.plot(pi_array_0_new,diff_new,label=r'high $\\&$ low $\\eta$')\n",
    "plt.plot(pi_array_0,diff2,label=r'full info $\\&$ low $\\eta$ ')\n",
    "plt.xlabel(r'$\\pi$', fontsize=16)\n",
    "plt.ylabel('$\\Delta \\Omega$', fontsize=16)\n",
    "plt.legend(fontsize=16)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "# DataOut = np.column_stack((pi_array_0_new[:-2],diff_new[:-2]))\n",
    "# np.savetxt('diff_highloweta_lowc.dat', DataOut)\n",
    "# DataOut = np.column_stack((pi_array_0,diff2))\n",
    "# np.savetxt('diff_fullinfoloweta_lowc.dat', DataOut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
